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I Statement of Objectives 


Our general objective for this work was to develop a theoretically and 
experimentally consistent explanation for the diffuse component of radar backscatter 
from Mars. The strength, variability, and wavelength independence of Mars 
diffuse backscatter are unique among our Moon and the terrestrial planets. This 
diffuse backscatter is generally attributed to wavelength-scale surface roughness and 
to rock clasts within the Martian regolith. Through the combination of theory 
experiment, we attempted to bound the range of surface characteristics that cou d 
produce the observed diffuse backscatter. Through these bounds we gained a 
limited capability for data inversion. Within this umbrella, specific objectives were. 


(a) To better define the statistical roughness parameters of Mars' surface so that 
they are consistent with observed radar backscatter data, and with the physical 
and chemical characteristics of Mars' surface as inferred from Mariner 9, the 
Viking probes, and Earth-based spectroscopy. 

(b) To better understand the partitioning between surface and volume scattering 
in the Mars regolith. 


(c) To develop computational models of Mars' radio emission that incorporate 
frequency dependent, surface and volume scattering. 


n Approach 

There is never enough radar scattering data to invert Mars radar backscatter 
profiles to obtain unambiguous estimates of surface roughness - even if we 
possessed the scattering theory that would allow us to do so. Our approach to t e 
problem has been to assume something about the statistical roughness 
characteristics of Mars' surface, and, based upon this roughness model, to develop a 
limited scattering theory or empirical law that could be used to relate observed 
backscatter to roughness. We are currently testing this approach by developing a 
classification algorithm based upon the empirical law developed in this 
investigation. Applying this algorithm to SIR-C data, we will attempt to classify 
Andean Volcanics according to their roughness. This application is funded by the 
SIR-C project and is not strictly a part of the Mars investigation. 


Ill Assumed Surface Statistics 

For roughness characteristics, we assumed a power-law surface. Several 
investigations have suggested that roughness profiles of many natural surfaces obey 
a power spectrum of the type. 




4 


Sz(f) = cf' P (1) 

[e.g., 7-8, 10, 22, 26, 28, 30-31, 34-35, 52, 60] where S z is the power spectral density, f is 
spatial frequency, c is a scaling constant, and P is a constant that is related to the 
fractal dimension of the surface [e.g., 21]. We tested this hypothesis on several flows 
of the Mount St. Helens volcano and found that fresh flows in that environment 
do, in fact, exhibit a power-law spectrum over several decades of spatial frequency 
both below and above the range of spatial frequencies of observing radars. In 
contrast, water deposited sediments of Mount St. Helens ejecta do not obey a power- 
law. While aeolian deposition on Mars may be like fluvial deposition at Mount St. 
Helens and not obey a power-law, the assumption that Mars volcanic deposition 
does obey a power law gives us a place to begin. We have presented various of these 
findings in several forums [65-66] and published the work in Geophysical Research 
Letters [2]. The GRL paper is included on the following 4 pages. 
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MULTI-SCALE ROUGHNESS SUIT I RA OK MOUNT ST. HELENS DEBRIS FLOWS 
Kirhani 1. Austin an-! Anthony W. Eneland 

Radiation Laboratory. Dopt. of El«*«*t i if.d i iif.'ri n si and Computer Science. University of Michigan 


\h>tnn!. A r« u : : : i » — -pei-trum allows surface struc- 
ture to he interprfied a- a <um of sinusoidal components 
with differtna wavrlrnat hs. knowledge of tin* roughness 
spectrum give*- insight into tlie meciianisnis responsible for 
elect romacnet ic scattering at a invert wavelength. Mea- 
sured sport ra from 1 0- vear-old primary debris flow sur- 
fares at Mount St. Helens conform to a power-law spectral 
model. -Mggesting that these surfaces are scaling over the 
measu, range of spatial frequencies. Measured spectra 
from water-deposited surfaces deviate from this model. 

Introduct ioti 

Recent investigations [Evans et al., 1988: Gaddis et ah. 
19S9] demonstrate that radar images are a useful tool in 
the study of volcanic terrain. Although there have been 
experimental investigations of radar scattering by such 
terrains, there have been few rigorous attempts to link 
the scattering cross section to quantitative surface rough- 
ness. Scattering from many volcanic terrains is difficult 
to model due to their extreme roughness and short cor- 
relation lengths. Rough surface scattering theories of- 
ten assume that surfaces possess a Gaussian autocorrela- 
tion function (and corresponding Gaussian roughness spec- 
trum) [Chen and Fung, 1988: Thorsos and Jackson, 1989]. 
In contrast, many natural surfaces are better described by 
power-law spectra [Berry and Hannay. 1978: Sayles and 
Thomas. 1978: Huang and Turcotte. 1 9S9, 1990: England, 
1992]. (These spectral forms are compared in Figure 1.) 
Rough surfaces that have structure over a wide range of 
spatial scales may he described using the concepts of frac- 
tal geometry introduced by Mandelbrot [1983]. Random 
rough fractal surfaces have power-law spectra: the slope of 
such a spectrum is related to the surface s fractal dimen- 
sion (Austin et al.. Special Problems in the Estimation of 
Power- Law Spectra as Applied to Topographical Model- 
ing. submitted to /ZTE7T Tvcmsactions on Ctcoscicncc und 
Rew off Sensing. 1993. hereinafter referred to as Austin et 
ah. submitted manuscript, 1993). In this report, we exam- 
ine estimates of the roughness spectra of volcanic debris 
flows at several spatial scales. 

Debris Flow Measurements 

Surface roughness measurements were performed on sev- 
eral debris flows near Mount St. Helens (Cascade Range, 
northwestern United States) during September 1990. Mt. 
St. Helens was chosen because of its ease of access and its 

Copyright 1993 by the American Geophysical Union. 
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extremely young terrain". I hi 1 >urta< e" that we examined 
were located in the debris avalanche weM • nort hwe>t ot the 
volcano, along the North fork loutle River valley i Fig- 
ure 2). Most of the debris was deposited during the erup- 
tion of Mount St. Helens on IS May l'N(). Since then, the 
deposits have undergone erosion by wind and water. 1 he 
sites chosen for measurement were visually representative 
of principal surface textures, accessible to logging roads, 
and located in large homogeneous areas which should be 
distinguishable in a radar image. Figure 3 shows one of 
the rougher surfaces. A geologic description of the debris 
avalanche is given by Glicken [1989], 

Our objectives were to characterize the surface rough- 
ness of the debris flows at scales smaller than, on the or- 
der of. and larger than the radar wavelength of common 
remote sensing radars. Two measurement techniques were 
used. A computer-driven, 2D laser profllometer recorded 
surface height profiles of square grids with sides between 
8 cm and 1 m in length. The grids were sampled at inter- 
vals ranging from 2 mm to 2 cm. We used surveying in- 
struments to measure elevation at grid points within larger 
square areas measuring 32 m on a side. Sampling intervals 
within the larger squares were 1-4 m. 

The 2D laser profllometer system was designed specif- 
ically for this experiment. Its main component is a sur- 
veying electronic distancemeter (EDM), which uses an in- 
frared laser to measure distance. The EDM is mounted on 
an XV table, allowing the EDM to scan a surface area of 
up to 1 m 2 . The EDM laser has a spot diameter of approx- 
imately 1.5 mm. In its most precise mode, the standard 
deviation of the measured surface height is 3 mm. Surface 
height profiles from a typical grid are shown in Figure 4. 
While the two large rocks are quite prominent, a larger 
sample of the surface would show even larger “lumps.” 
Larger-scale topography was surveyed. Square grids 
with sides of 32 m were delineated by cables with markers 
at 1 m intervals. A self-leveling surveying level and sta- 
dia rod were then used to measure the surface height at 
intervals of 1-4 m, depending on the surface roughness. 
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[ i<l. 4. The .Johnston Ridge site (JRS) is located about 
10 km northwest of the crater. I he Elk Rock site (ERS) 
is located farther west, about 10 km west -northwest of the 
crater. 


Spectral Estimation 

Estimates of the roughness spectrum were calculated 
from both types of elevation data. Profilometer data suf- 
fered from errors which prevented the direct application 
of spectral estimation techniques. Since the measurements 
could not be repeated, it was necessary to apply correct i\e 
filtering. The procedure is fully described in Austin and 
England [1991]. 

As previouslv noted, investigators have found that man\ 
natural surfaces are scaling. Linear profiles of such surfaces 
have (over the range of measurable spatial frequencies) 
a surface roughness spectrum (power spectral density of 
surface height) >z that has the form of a power law: 

$ 7 .( f ) — c l/l ^ 


where / is the spatial frequency and c and j3 are constants, 
with 1 < 3 < 3. 3 is sometimes called the spectral slope; a 
power-law spectrum is linear with slope ( — 3) when plotted 
on a log- log scale. 

Power-law surfaces with isotropic statistics have a two- 
dimensional roughness spectrum Sz2D that has a similar 


form: 

$ Z 2 D { fr ) = a fr ^ 


( 2 ) 


where f r is the radial spatial frequency satisfying f r 
p + p 9 and 2 < -y < 4. Because estimation of direc- 
tional spectra requires a much greater quantity of data, 



Fig. 3. A sample surface scanned by the EDM. The sur- 
face lias been spray- painted to increase its reflectivity. 
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we assume in this study that the surfaces have isotropic 
statistics as a first -order approximation. 

Power-law spectra introduce unique difficulties in the 
spectral estimation process. In a related paper (Austin et 
al.. submitted manuscript. 1993). we describe how ( apon s 
estimator provides an estimate of the roughness spectrum 
that is better than that given by Fourier-based estimators 
(such as the periodograrn ) in the power-law case. ( apon s 
estimator (described in Kay [1988]) essentially customizes 
a filter at each frequency of interest to minimize the total 
power output subject to the constraint that the gain at 
the frequency of interest is unity. In the power-law case, 
the low-frequency side lobes are reduced, thus Capon s es- 
timator suffers less from spectral leakage, which can make 
Fourier estimates insensitive to the power-law slope 3. 
Capon's estimator also has less variance than comparable 
Fourier- based estimators. 

Capon's estimator, Pcap(I)' * s obtained by first cal- 
culating an estimate of the autocorrelation matrix R zz< 
whose elements are defined as 


[R zz]tci — E [Z[n]Z[n + k — /]] , (3) 


where Z[n] is the surface height and £[•] indicates an ex- 
pected value. Capon’s estimator is then given by 


Pc.APif) 


PA 

e w Rzz e 


(4) 


where p is the maximum dimension of the autocorrelation 
matrix, A is the sampling interval. Rzz i nverse 

an estimate of the autocorrelation matrix, and 


| c j2 rr/A e j4r/A 


e }2rr{p-\)J* j 7 ^ (5) 


in which / is the spatial frequency. e H is the Hermitian 
transpose of (5). 

The bias of Capon’s estimator is independent of N (the 
number of points in a sampled row) but is dependent on 
p, the dimension of the autocorrelation matrix. For a data 
set of length A r . an increase in p will result in reduced bias 
at the cost of increased variance. We set p ~ 0.3 A and cal- 
culated a Capon's spectral estimate for each row of each 
profilometer or survey grid. (The mean row height was 
subtracted before calculating the estimate.) The spectral 
estimates were t hen averaged over the rows of each grid. 
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Fig. 5. (A) Survey and profilometer spectra from a JRS 

primary debris flow. (B) The same survey spectrum with 
a profilometer spectrum from an adjacent debris flow. (C) 
Spectra from a mixed terrain survey and a primary debris 
flow scan at ERS. 


The high-frequency half of each averaged estimate was dis- 
carded to reduce the effect of aliasing. Averaged estimates 
at a single scale or at multiple scales were then fit with a 
power-law function (1), where appropriate, using a min- 
imum absolute deviation fit. Profile (ID) spectra fitting 
the power-law model were then converted to surface (2D) 
spectra using the formulation given in (Austin et ah, sub- 
mitted manuscript. 1003). in which the ID and 2D rough- 
ness spectra are related through their corresponding auto- 
correlation functions and the 2D autocorrelation function 
is assumed isotropic. 

Although this method of determining the surface spec- 
trum from linear profiles was made necessary by the pro- 
filometer problems, it is not without benefit. One grid of 
measured surface heights yields many linear profiles. The 
profile rows are not completely independent, but the row 
spectral estimates can still be averaged to obtain a signifi- 
cant decrease in variance, allowing us to set p (in Capon's 
estimator) higher to decrease the bias. This averaging 
would not be possible with a two-dimensional estimator, 
given the same measured data. 
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of material hv wind and water. Spectra tn>m the romMi 
est of t he-e are shown in Figure a A. I he survey -peel mm 
is the curve extending from 0.(1’ to U.23 m“ l ; the higher 
frequenrv spectra were computed from profilometer scans. 
The linear character of these spectral estimates suggests 
that a power-law spectrum is appropriate tor this surface 
The equation of the power-law lit to these estimates is 

>■;,(/) = (3.21 X 10' 4 t;/r J H V (fi) 

(The hat denotes an estimate.) This fit is shown as a 
dashed line in f ignre oA. W hile we do not know that the 
spectrum follows t his power law at the intermediate (un- 
measured) spatial frequencies, we do think that il is sig- 
nificant that both the meter- and centimeter-scale spectra 
arc well fit by the same power law. 

We can convert (6) to an estimate of the surface (2D) 
spectrum, obtaining the following power law: 

$Z 2 o(fr) = (E77 x 10 4 )/ r • (*) 

In Figure 5B, the same survey spectrum is shown to- 
gether with a profilometer spectrum collected on a similar 
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Fig. 6. (A) A profilometer spectrum (A = 1 cm) of an 

ERS sedimented surface, shown with the survey spec- 
trum of the local mixed terrain, the power-law fit (10), 
and the EDM noise floor. (B) A profilometer spectrum 
(A = 5 mm) of a .IRS sedimented surface, shown with the 
survey spectrum of the surrounding primary debris flow, 
the power-law fit (fi). and the EDM noise floor. 
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fit mav hr appropriate here as well and is shown in the 
figure*. I hr equations lor this fit are 

S>(/> = ( 4.9$ x 10- 4 )!/r 2M (10) 

Szmi.fr > = (2.86 x 10 - * 1 )/~ 3 51 (11) 

Surfaces formed of water-deposited sediments within the 
hummock v debris had spectra that were markedly differ- 
ent from those of the primary surfaces. For example, in 
Figure 6 A. we see that the spectral estimate from a scan 
of such a surface lies on the noise floor due to the lim- 
ited precision of the EDM. The true spectrum of this sur- 
face must lie below this level and therefore cannot form 
a power-law trend with the survey spectrum of the local 
mixed terrain. This difference in spectral behavior is con- 
sistent with the visible character of the sedimented sur- 
face, which appeared quite flat. Similar remarks apply to 
the profilometer spectrum of a water-deposited sediment 
at the .Johnston Ridge site, which is shown in Figure 6B 
with the survey spectrum of the surrounding primary de- 
bris flow. 

Conclusions 

Spectral estimates derived from multi-scale topographi- 
cal measurements of several debris flow surfaces at Mount 
St. Helens show differences that correspond to observed 
surface characteristics. The extremely rough surfaces affil- 
iated with 10-year-old primary debris flows possess rough- 
ness spectra that are well modeled by a power-law func- 
tion at certain length scales, implying that these surfaces 
are scaling over the measured range of spatial frequencies. 
Surfaces created through sedimentation have spectra that 
lie below the noise floor at higher spatial frequencies and 
tints do not fit the power-law model. 

While we are not vet able to develop a classifier of vol- 
canic terrains based on surface roughness spectra, this 
studv shows that a power-law spectrum is a reasonable 
descriptor of some volcanic surfaces. The scattering prop- 
erties of such surfaces are not well understood. There has 
been some application of standard scattering models to 
surfaces of this type, but rougher volcanic terrains often 
violate assumptions made in these models (e.g., rms height 
< radar wavelength). The present study and the others 
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IV Estimation of Power-Law Spectra 

Having measured topography at scales from millimeters to 32 meters on 
several volcanic flows, our next task was to estimate the roughness spectra of these 
topographies. We immediately found that simple Fourier-based estimators, like the 
periodogram, produced greatly distorted spectra. The low frequency leakage of a 
power-law surface produced an artificial [3 in equation (1) that tended to lay near 2. 
Many of the early papers on the spectra of natural surfaces [e.g., 1] have suggested 
that the p for natural surfaces always lay near 2. We believe that this conclusion 
may have been the result of the leakage phenomenon that we observed. In fact, 

while the Mount St. Helens debris flows obeyed a power-lay, their p always lay 
between 2.31 and 2.51 - which are significantly different from 2. An artificial surface 

produced with a p of 2 would appear too rough at the higher spatial frequencies 
when compared to the Mount St. Helen surface profiles. 

We presented these findings at the Lunar Science Conference [67] and 
published the results in the Transactions of Geoscience and Remote Sensing [3]. The 
TGARS paper is included on the following 12 pages. 




iF.fcf- TRANSACTIONS OS UH’St. II SO \Mi RTMoTI \ Ol 


V ■ i 


HO. 


Special Problems in the Estimation of Power-Law 
Spectra as Applied to Topographical Modeling 

Richard T. Austin. Student Member , IEEE, Anthony W. England. Senior Member, IEEE, 
and Gregory' H. Wakefield, Member, IEEE 


Abstract — An increasing number of topographical studies had 
that natural surfaces possess power-law roughness spectra. 
Power-law spectra introduce unique difficulties in the spectral 
estimation process. We describe how an improper window 
choice allows leakage that yields a spectral estimate that is in- 
sensitive to the spectral slope. In addition, the commonly used 
Fourier-based spectral estimates have higher variances than 
other available estimators. Higher variance is particularly 
problematic when data records are short, as is often the case 
in remote-sensing studies. We show that Capon’s spectral es- 
timator has less variance than Fourier-based estimators and 
measures the spectral slope more accurately. We also show how 
estimates of a 2-D roughness spectrum can be obtained from 
estimates of the 1-D spectrum for the isotropic power-law case. 

I. Introduction 

S PECTRAL analysis is a tool of increasing importance 
in the characterization of natural surfaces. One- and 
two-dimensional spectral estimates are useful because they 
allow surface structure to be interpreted as a sum of si- 
nusoidal components with differing wavelengths. A num- 
ber of studies have found that spectra computed from one- 
dimensional surface profiles may be reasonably modeled 
using a power-law spectrum of the form 

S z (f) = c\f\~ 0 (1) 

over some range of spatial frequencies. S z (f) is the power 
spectral density of the surface height random process; it 
has units of meters 2 /meters“ 1 or meters 2 per unit spatial 
frequency. The spectral exponent 0 indicates the relative 
contribution of different wavelength components. Higher 
values of 0 indicate that the surface roughness is mostly 
due to low-frequency components, while lower values de- 
note significant high-frequency contributions. 0 is also 
known as the spectral slope; a power-law spectrum is lin- 
ear with slope — 0 when plotted on a log-log scale. The 
roughness amplitude c is a multiplicative factor scaling 
the roughness at all spatial frequencies. Applications of 
power-law spectra include studies of regional topography 
[l]-[5], rock surfaces [6], seafloor morphology [7]-[ll], 
and other surfaces [12], [13]. 

Rough surfaces that have structure over a wide range 
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of spatial scales may also be described using the concepts 
of fractal geometry introduced by Mandelbrot [14]. Frac- 
tal objects are continuous but not differentiable; the frac- 
tal dimension D f is a real-valued measure of how a line 
(surface) fills a plane (space). For a one-dimensional 
(1-D) surface profile. D f takes on values between 1.0 
(smooth and differentiable) and 2.0 (plane-filling). Two- 
dimensional (2-D) surfaces have D f between 2.0 and 3.0. 
Adler [15] shows that the surface profile created by the 
intersection of a plane and a 2-D fractal surface is itself 
fractal with a fractal dimension equal to that of the 2-D 
surface decreased by one [2], Random rough fractal sur- 
faces have power-law spectra; Mandelbrot and Van Ness 
[16] and Voss [17] derive the relation between and the 
spectral slope 0 of a 1-D profile of a rough fractal surface: 


Thus, the spectral exponent satisfies 3 > 0 > 1 for 1 < 
D f < 2. 

Although synthetic topographies appear most realistic 
with D f * 2.2 [14], [17], yielding linear profiles of D f = 
1.2 and 0 = 2.6, several studies [1], [4], [7], [8] have 
reported measured values of $ (the hat denotes an esti- 
mate) near 2.0, leading some investigators to suggest that 
spectral slopes of 2.0 are characteristic of natural topog- 
raphy [1]. One explanation for the discrepancy between 
measured £ and inferred values of 0 is that spectral anal- 
ysis of nonstationary regions tends to result in 3 near 2.0 
[9], [18]. We show in this paper that leakage in spectral 
estimators may also yield $ that cluster near 2.0. 

While leakage is a well-known problem of classical 
spectral estimators such as the periodogram, common 
Fourier-based spectral estimators also suffer from high 
variance. The usual remedy is to average spectral esti- 
mates derived from multiple data sets or multiple seg- 
ments of a long data set. This is not always possible: in 
remote-sensing studies, surface profiles are often only one 
of many “ground- truth” parameters to be collected. As a 
result, investigators may be faced with the problem of es- 
timating the surface spectrum using a limited data set. 1 In 

'In the case of 1-D surface profiles, we find that a typical record consists 
of at most 500 points; 2-D profiles are even shorter, rarely exceeding 50 
x 50 points. Numbers of this magnitude are considered small from a sta- 
tistical inference viewpoint, where samples of thousands or hundreds of 
thousands of points are more common. 

00 © 1994 IEEE 
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this paper, we describe a spectral estimator that reduces 
the leakage problem and predicts c and 3 with reduced 
vanance in the case of short data runs 


II. Spectral Estimates from Linear Profiles 

Fourier-based spectral estimators le g., the periodo- 
aram) render an estimate of the power spectral density 
whose expectation or mean value is a convolution of the 
spectrum of the sampled profile and the Fourier transform 
of a window function introduced by the finite extent of 
the sample profile. Various window functions are used: 
common windows include rectangular, triangular, and 
Hanning windows. The selection of a particular window 
is based on a trade-off between spectral resolution and 
spectral leakage. 

Spectral leakage refers to the inaccuracy of a spectral 
estimate at a given frequency due to convolution with the 
transformed window; i.e., spectral power “leaks’ from 
nearby frequencies according to the shape of the trans- 
formed window. We show spectral domain window func- 
tions corresponding to a periodogram and a modified pe- 
riodogram with Hanning window in Fig. 1 . At each point 
or frequency at which the spectrum is to be estimated, the 
true spectrum (the dashed line in Fig. 1) is weighted by 
the window function centered on that frequency and then 
averaged over frequency. Contributions from frequencies 
other than the frequency of interest constitute spectral 
leakage. We can reduce spectral leakage by choosing an 
estimator whose window has sidelobes that decay quickly 
(such as W H in Fig. 1), but such a window will have a 
broader main lobe, decreasing the spectral resolution. (See 
[19] for a discussion.) In the power-law case, leakage 
from the large peak at low spatial frequencies affects the 
estimate at upper spatial frequencies (making the slope 
shallower) unless a window function with very low side- 
lobe levels is selected. We now show this result for the 
case of 1-D spectral estimates based on 1-D surface pro- 
files. 

Spectral estimates derived from sampled data suffer 
from aliasing if the sampled process has spectral com- 
ponents at frequencies greater than the Nyquist frequency 

/o 



where A is the sampling interval. In the sections that fol- 
low, we will assume that aliasing is not present. Simu- 
lated data sets will be bandlimited to prevent aliasing. 
Modifications due to aliasing will be presented in a later 
section. 

A. Fourier-Based Estimators 

Consider a one-dimensional surface Z(x), where Z, the 
surface height, is a single-valued function (i.e., there are 
no overhangs). Suppose that the surface height has been 
measured at N locations x, (/ = 0. • • • , N — 1), which 



Fig 1 Spectral domain window functions M ( l(h tor the periodogram 
and W H (13) tor the modified periodogram with Hanning window, both 
with V = 25b. The convolution ot the window Junction with the true spec- 
trum le g . the power-law spectrum shown as the dashed line) ma> result 
in a spectral bias lan inaccuracv in the mean value ot the estimate > 


are spaced at intervals A. The periodogram spectral esti- 
mator P PER (/d f° r the sampled surface height is given by 
Kay [20]: 


P PER(/a) “ yy 


,v - i 

Z Z[.r„] exp ( ~j2irf±n) 

n - 0 


( 4 ) 


The frequency is written as a reminder that (4) is given 
in terms of a normalized frequency/^, where 

h = ^ ( 5 ) 

and f A satisfies -y < / A < y. (Most spectral estimation 
texts use normalized frequency; we will follow this con- 
vention until a later section when we compare spectra ob- 
tained from sampled profiles with differing A.) 

While the periodogram may be evaluated for any |/J 
less than y (the normalized Nyquist frequency), it is often 
calculated using a fast Fourier transform (FFT) and con- 
sequently is evaluated only at discrete frequencies = 
i/N, where / = —N/2, ■ • • , -1,0, 1, ■ • • , Nil. To 
evaluate the periodogram using an FFT at additional fre- 
quencies, we can “zero-pad” the original data senes to 
length M by adding M - N zeros to the end. The period- 
ogram will then be evaluated at frequencies = HM. 
Most FFT algorithms require that N be an integer power 
of two. Series for which N is not an integer power of two 
must be zero-padded up to the next higher power of two. 

The periodogram is an unreliable estimator of the power 
spectral density because it has a variance that is equal to 
the square of its expected value, independent of N [20]. 
The usual practice in the field of spectral estimation is to 
average multiple periodograms to reduce the variance, but 
this may not be possible if the quantity of data is limited. 2 
We show later in this section that P PER is a poor estimator 
in the power-law case even if the variance problem were 
absent. 


2 For example, suppose the data record is 1000 points long. To reduce 
the variance by a factor of eight, we can divide the data record into eight 
segments of length 125, calculate eight spectral estimates, and average the 
results. The bias of the estimate will be increased (due to the broader main 
and side lobes of the window), and the reduction in variance may be less 
than eightfold if the segments are not uncorrelated. If the entire data record 
has only 125 points, this procedure is not useful 
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A modified penodogram /> HAN ( A> may be obtained by 
multiplying the data senes by a Hanning window before 
calculating the spectral estimate: 

1 | v -‘ r 

>WA> = 777 ;! — Z[x n ]\v H M exp i -j2irf A n) (6) 

A L \ n = II 

where the Hanning w indow 

1 i llvtl \ j 

"h\”\ = ^ 1 - ‘■'os y J j ,7) 

replaces the rectangular window’ ( v\ > [ n] = 1) implicit in 
(4), and an extra normalization factor is introduced: 

1 V_1 

u = T, £ Wtf[« 1- (8) 

N n = 0 

(Equation (6) follows the formulation of Welch [21] using 
a single data segment.) The transform of the Hanning 
window w H [n] has sidelobes that are lower and a main 
lobe that is roughly twice as wide as that of the transform 
of a rectangular window w R [n] of the same width. Con- 
volution with a transformed Hanning window results in 
less leakage from low spatial frequencies (because the 
sidelobes are smaller) at the cost of decreased spectral res- 
olution (smoothing due to the wider main lobe). We must 
consider the smoothness of the spectrum and how quickly 
it rolls off and then decide whether the leakage reduction 
is worth the loss in resolution. In the power-law case, 
reducing the leakage seems more important. 

We now examine the expected values of these two es- 
timators for a power-law surface spectrum. The expecta- 
tion of the periodogram is given by [20] 

J l/2 

(9) 

-1/2 



normalized spatiaJ frequency, dimensionless 


Fig . expected value of the penodogram spectral estimator for power 
law sunaces with J = 2.0 and 2.4 and t ^ =1.0- 10 f '. The estimates are 
based on 256-point profiles that were zero-padded to 16 384 points so that 
the penodogram would be evaluated at mans frequencies. Solid dots indi- 
cate the expected values of the penodogram when evaluated without zero- 
padding. 


where 


W h (L) = — — (b, + - b, + - b x 
H 4 NU V 2 - 2 " 


bi 


sin (x/ A /V) 
sin (x/ A ) 


sin x 


b , = 


h - 


N - 1 


N 


sin x 


sin < x 


b , = 


L - 

/a + 


1 


N - 1 
1 

N - 1 


N 


sin i x 


/a + 


N - 1 


(13) 

(14) 


(15) 


(16) 


where W^/a) is the transform of a triangular window: 


WtU. A) 


1 / sin x/ A *y 
N \ sin x/ A / 


( 10 ) 


and S ZA (/ A ) is the normalized form of the true power spec- 
tral density of Z(x) chosen such that J Sza(/a) df A = 
\SAf)df- 


5za(/A) = | SAf) 



= A 0 1 c(fA) 8 

5za(/a) = c A fl 0 . (11) 


The expected value of the modified periodogram esti- 
mator P HAN has a form similar to (9) [23, p. 553]: 


han(A)1 



W h (/a ~ t)SzA(Z)dZ 


( 12 ) 


and U is given by (8). 

To show the performance of these Fourier-based esti- 
mators, we calculate their expected values using (9) and 
(12). We avoid the singularity at the low end of a pure 
power-law spectrum by modeling the exact surface spec- 
trum as a Rayleigh function at low spatial frequencies: 


f (ala 2 ) | / A | exp(-/i/[2a 2 ]), 


Sz a(/a) = { 


I / A I SS 0.0001 

c aI/aI B 

0.0001 < |/ A | < J 


(o, I/aI > j 


(17) 


where the parameters a and a are chosen by enforcing the 
continuity of S ZA and its first derivative at |/ A | = 0.0001 . 

Fig. 2 shows the expected value of the periodogram 
spectral estimator calculated from theoretical spectra with 
spectral exponents of 2.0 and 2.4. The oscillatory behav- 
ior is due to the sidelobes of W T (f±). These oscillations 
are not visible if the periodogram is evaluated only at fre- 
quencies / = UN (i.e., without zero-padding) as indicated 






Fig 3 Expected value of P PER for five 256-point power-law profiles eval- 
uated at frequencies corresponding to the case ot no zero-padding. 


in the figure. Since it is difficult to fit power-law functions 
to such oscillatory estimates, we sampled the expectations 
of periodograms of surfaces with five different spectral 
exponents at frequencies corresponding to the case of no 
zero-padding (Figure 3) . 3 These curves may be compared 
to the exact spectra in Figure 4. 4 While the exact spectra 
vary in slope, the periodogram estimates are insensitive 
to the slope for several values of /3. Fitting power-law 
functions to the linear portions of these estimate expec- 
tations (f A < 0.2), we see (Table I) that the estimated 
spectral exponent ; clusters near 2.0 for 2 < /3 < 3. 
These values of 13 correspond to fractal dimensions that 
are typical for natural terrains (1.0-1. 5 for profiles, 2.0- 
2.5 for surfaces). 

Slope insensitivity is not a problem with Phan(/a)> 
whose expected value is shown for various /3 in Fig. 5. 
The expected values of these estimates closely approxi- 
mate the exact values for frequencies at which the main 

3 We emphasize that Fig. 3 shows expected values, £[^per(/)]. 0-C., 
mean values) of the periodogram estimator, calculated using (9) at fre- 
quencies corresponding to the case of no zero-padding for random profiles 
Z(x) having power-law spectra given by (17) with five different values of 
/3. Plots of the estimator P PER itself would show wide variance and depar- 
tures from linearity (similar to the estimates shown in Fig. 8) and would 
vary for different realizations of the random process, i.e., for different sam- 
ples of a given surface. 

The expected values of the estimator appear linear at the low-frequency 
end because we have assumed in (17) that the true surface spectrum is 
power-law (with linear slope) down to a frequency of/ A = 0.0001 — a fre- 
quency too low to be accurately measured by the sampled profile segment. 
We used this assumption because surface profiles are often too short to 
show very long wavelength roughness components. (This was certainly true 
for our field work.) 

4 The expected values of the periodogram estimator in Fig. 3 have a higher 
total energy level than the corresponding exact spectra in Fig. 4 due to the 
spectral leakage that is also responsible for the slope insensitivity. Remem- 
ber that £[P per ] is the convolution of the true spectrum S ZA with the trans- 
form of a window function introduced by the finite surface sample [W T in 
(9)]. If the estimator were perfect, W T would be a delta function, and the 
estimate at frequency / would reflect only power in the true spectrum at 
that frequency. In the periodogram case, however, W T is given by (10) (a 
sine 2 function), so power from many frequencies above and below /affect 
P PER (/). Leakage from lower frequencies is particularly significant in the 
power-law case. 

5 For a given single periodogram, |3 may indeed take values between 2.0 
and 3.0, due to the variance of the periodogram. Long data records do not 
reduce this variance; the variance is independent of N (see [20, p. 422]). 
One may reduce the variance by segmenting the long data record as de- 
scribed previously, but this process will lead to the expected values shown 
in Fig. 3. Sampling P PER at frequencies other than f Al = i!N could lead to 
a variety of incorrect slopes, as can be seen in Fig. 2. 



Fie. 4. Exact values of the power-law spectra used in calculations of the 
expected values of the spectral estimates. 


TABLE I 

Power-Law Parameters Derived from Fits to the Expectations of 
Periodogram (1 A < / A < 0.2). Modified Periodogram <4 A r < / A 
< 0.5). Prewhitened H >N < / A < 0.2), and Capon s ( 1/(2/?) < / A 
< 0.5) Estimators 


Theoretic*! 

e[Ppeh] 

E[Phas] 


E[Pcap\ 

c A = t o to-* - 
0 = 1.2 

e A =9 71 10"' 
0- 1.231 

c A = 9 99 10"' 
0 = 1.201 

c A = 8 42 10“ ' 
0 = 1 255 

c A - 1 01 
0 = 1 198 

c A = l 0 10“° 
0 = 16 

c A = 1.04 10” 
0 = 1.656 

c A = 9.99 10- T “ 
0 - 1 601 

c A = 8 63 10" ' 
0 = 1.643 

c A = 1.01 10- tt 
0 = 1.598 

e A = 1.0 10“° 
0 = 2.0 

c A = 2.10 10‘ fl 
0 = 1.982 

c A = 9.96 ltP 
0 = 2.001 

c A = 8 66 10" ' 
6 = 2.042 

c A = 1.02 lO" -5- 
0 = 1 998 

c A = 1.0-10- 
3 - 2.4 

c A = 1.87 10- t_ 
0 = 2.032 

e A = 9.98 10"' 
0 = 2.402 

c A =8.65 10“' 
0 = 2 443 

c A = 1.02 10-"* 
0 = 2.397 

c A = 1,0 10'° 
0 = 2.8 

c A =3.83 10" 
9 = 1.986 

c A = 9.96 10- ’ 
0 = 2.804 

c A = 8.43 fiF 7 
0 = 2.859 

c A = 1 .02 10 _ * 

0 = 2 796 


lobe and first two sidelobes do not overlap the low-fre- 
quency peak (/^ > 0.015625). (Power-law fits to the up- 
per regions of these estimates are listed in Table I). While 
P HAN has a mean value that closely approximates the exact 
spectrum, this estimator suffers from a deficiency that is 
common to Fourier-based estimators: undesirably high 
variance. We examine the variance of P han(/a) * n a later 
section. 6 

B. Prewhitening 

The periodogram (Fig. 2) and, to a lesser degree, the 
modified periodogram (Fig. 5) suffer from spectral leak- 
age in the power-law case, especially for spectra with 
2 < 0 < 3. Fox and Hayes [9] and Gilbert and Malin- 
vemo [10] avoid the leakage problem by using pre whiten- 
ing procedures in which they modify the surface height 

''We have also examined the Blackman-Tukey spectral estimator P^, as 
described in Blackman and Tukey [22]. For processes other than white 
noise, the BT estimator decreases variance but increases bias as compared 
lo the periodogram (Kay [20, p. 80]). We did not include Porr in our com- 
parison because it is even more susceptible to spectral leakage (and there- 
fore, slope insensitivity) than P PER . The expected value of P B t is equal to 
the expected value of P PER convolved with another window, as derived in 
Kay [20, p. 98]: 

S l/2 

mh - o^i/ws)] ft 

- 1/2 

Since the BT estimator has an expected value that is the result of the con- 
volution of the true spectrum with two window functions, its value at a 
given frequency can be greatly corrupted by leakage from lower frequen- 
cies. 
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Fig. 5. Expected values ol P HAS ( /^i for five different 256-point pouer- 
law profiles. The profiles were zero-padded to 16 384 points. 



Fig 6 Expected values ot \P pu { f±) for hve different 257-point power-law 
profiles. The profiles were zero-padded to 16 384 points 


dan to have a relatively flat spectrum, perform spectral 
estimation, and then correct the estimate according to the 
original modification. 

We will examine here the simplest prewhitening ap- 
proach similar to that used by Gilbert and Malinvemo. If 
the power spectral density of the surface height is defined 
as [24, p. 270] 

S z (f)= lim^fHf Z(x) exp (-j2wjx) dx 
7“- oo / LI J -772 


(18) 

then the derivative of the surface height dZ(x)/dx will have 
a related spectrum: 

$z'(f) = lim -e\\ f J ^ exp (-y 27r/r) dx 

T-o>T Li -772 dx 

i n f 772 2_ 

- lim -E\ I/2t/ I Z(x) exp (—jlicfx) dx 

T -* oo I LI J — 77 2 

= 4ir 2 / 2 S 2 (/) 

= 4tt 2 c\ f\~ a + 1 

S z ^f) = c'\f\- 3 ‘ (19) 

and we see that the spectral exponent /3 ' of the derivative 
process takes on values between - 1 and 1 for 1 < 0 < 
3. The spectrum of the derivative process is more nearly 
flat (it is white noise for 0' =0). Since there is no pro- 
nounced peak at low frequencies, spectral leakage is less 
of a problem. 

Gilbert and Malinvemo approximate the derivative by 
taking first differences of the sampled data: Z'[x,] * 
Z[x, + ,] — Z[Xj], They then apply a Hanning window and 
use a periodogram to obtain a spectral estimate of S z - A (/ A ). 

After averaging multiple periodograms, an estimate of the 
spectral exponent # is obtained from an estimate of $ 

& = & ' + 2. (20) 

We now examine this procedure for the limited data case 
in which multiple periodograms are not available for av- 
eraging. We omii the Hanning window for simplicity. 

We define a prewhitened spectral estimate P PW (A) in 


terms of the first differences of the sampled data: 


^pw(A) — — 


.v - i 


^ {Zl x n+i] ~ Z[x n ]} exp i-j2irf A n) 


(21) 

The expected value of the prewhitened estimator may be 
derived by a procedure similar to the derivation of (9): 

pw(A)l = 2 ( W/C/a — 0 $ za (£) 

J — 1/2 

• [1 - cos (2 tt£)] (22) 

Using (22) and the modified power-law spectrum defined 
in (17), the expected values of the prewhitened estimator 
are calculated and are shown in Fig. 6. These estimates 
are based on 257-point profiles. 

Power-law functions were fit to the linear portions of 
these estimate expectations (/ A < 0.2), and the parame- 
ters c' and were converted to c and # using (20) and 
c = c'/(4ir 2 ). The results, given in Table I, show that the 
expected value of the spectral slope is fairly accurate be- 
cause the spectral leakage has been largely eliminated. A 
Hanning window may still be useful in practice (since we 
do not know a priori that a spectrum is power law), but 
we see that leakage is not a problem for the prewhitened 
estimator in the power-law case. The variance of P pw (/ A ) 
will be examined in a later section. 


C. Capon ’s Estimator 

An alternative procedure that is more direct than the 
prewhitened periodogram is the use of Capon’s estimator 
[25] (the so-called “minimum variance spectral esti- 
mator”) described in Kay [20, chap. 11]. This estimator 
was developed originally for geological signal processing 
problems in which the number of sensors (and therefore 
the number of spatial samples) were extremely limited and 
posed restrictions on what could be inferred from standard 
spectral analysis. Capon’s estimator essentially custom- 
izes a filter at each frequency of interest to minimize the 
total power output, subject to the constraint that the gain 
at the frequency of interest is unity. Therefore, the filter 
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may be asymmetric in sidelobe level according to the 
shape of the signal spectrum; in the power-law case, the 
sidelobes are adjusted to reduce the leakage from low- 
frequency components. 

Capon’s estimator P c ap< * s obtained by first calcu- 
lating an estimate of the autocorrelation matrix R^, whose 
elements are defined as 

[Rzz], = E\Z*\n)Z[n + i - j}\. (23) 


We can estimate the elements of the autocorrelation ma- 
trix using the modified covariance method, described in 
1201 : 


1 


(Rzzl-7 


2 (TV - p ) 


Z Z\n - i]Z[n - j] 


rV - 1 - p 

+ Z Z[n + i]Z[n 

n = 0 



Capon's estimator is then given by 


^cap(/a) 



(24) 


(25) 


where p is the dimension of the autocorrelation matrix and 
e = [1 e ilrfx e ,4rh ••• e l2 * p ~ Wi ) T . (26) 

Since (24) is an unbiased 7 and consistent estimator of 
the covariance matrix, we can evaluate the expected value 
of Capon’s estimator by substituting the exact covariance 
matrix (obtained from the exact theoretical spectrum) into 
(25). This procedure was used to calculate the expected 
values of P CAP (f A ) for five values of 0 using a covariance 
matrix of dimension p = 70. The expected values of these 
estimates (Fig. 7) have good agreement with the exact 
spectra over a wide range of spectral slopes. Parameters 
of power-law functions fit to the expected values of 
P CAP (f A ) for f A > l/(2p) are compared with the exact 
spectral parameters in Table I. 

The bias of Capon’s estimator is independent of N be- 
cause (24) is an unbiased estimator of ftzz- However, the 
bias does depend on p ( pA is the longest lag for which 
the covariance is estimated). For a given data set, an in- 
crease in p will result in reduced bias at the cost of in- 
creased variance. Thus, N indirectly influences the bias 
because the range of p is constrained by N. We generally 
found p ~ 0.3 N to be a good compromise between vari- 
ance and bias. 

In the power-law case, the bias of Capon’s estimator 
was found to increase noticeably for spatial frequencies 
corresponding to wavelengths much longer than pA. We 
therefore discard calculated values of P cap (/a) f° r spatial 
frequencies below l/(2p). 


’Equation (24) is approximately but not strictly unbiased because all lags 
are not weighted equally and overlapping lags are used. A variety of co- 
variance estimators have been examined, and (24) was found to be the 
preferred estimator for short data records. See Kay [20] or other texts for 
more information. 



Fig 7. Expected values of P t » p l /»> for synthetic profiles with various 3 
usine an autocorrelation matrix of dimension p = 70 Estimates are eval- 
uated at the same 16 384 frequencies used in Figs 2. 5. and 6 


D. Variance Comparison 

Capon’s estimator, the modified periodogram with 
Hanning window, and the prewhitened periodogram all 
produce parameter estimates whose expected values 
closely approximate the exact values for a power-law 
spectrum. We now compare the estimators in terms of 
variance. The exact expression for the variance of the 
modified periodogram involves higher order moments and 
is usually evaluated only for a Gaussian white-noise case 
[21]. The statistical properties of Capon’s estimator for 
time series data are not known [20] , although some results 
have been derived for array data. To avoid these problems 
and test these estimators under actual-use conditions, we 
generated short synthetic topographic profiles of known 
statistics using a spectral synthesis algorithm. We then 
compare parameters obtained from Capon’s, the modified 
periodogram, and the prewhitened periodogram estimates 
with those of the spectra used to create the profiles. 

Our spectral synthesis algorithm closely follows [26] 
and consists of the following steps. 1) Generate a set of 
discrete Fourier amplitudes that, when squared and mul- 
tiplied by l IN, satisfy the desired power law. While [26] 
used a pure power-law spectrum, we use the modified 
spectrum (17) to assure a zero-mean surface height. 2) 
Multiply the amplitudes by a Gaussian random variable 
such that the mean value of the amplitude satisfies the 
power law. 3) Generate a set of complex discrete Fourier 
coefficients (i.e., the FFT of a real surface) using the ran- 
domized amplitudes and a uniformly distributed random 
phase, enforcing symmetry conditions such that the in- 
verse FFT will be real-valued. 4) Calculate the inverse 
FFT of the coefficients, resulting in a synthetic surface 
profile. 

Using this algorithm, we generated 10 synthetic 64- 
point surface profiles with 0 = 2.4 and c A , a, and o cho- 
sen as before. We then calculated 1) Capon’s estimates 
using a covariance matrix of dimension p = 20, 2) mod- 
ified periodogram estimates using a full-width Hanning 
window, and 3) periodogram estimates based on the first 
differences of the surface profiles. We see in Fig. 8 that 
the modified periodogram and prewhitened periodogram 
estimates have noticeably greater variance. 
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Fig. 8 Capon s, modified periodogram, and prewhuened periodogram 
spectral estimates tor 10 different power-lav* profiles The profiles were 64 
points long The modified periodogram and Capon's estimates have been 
offset tor clants The solid lines represent the exact spectrum. 


TABLE II 

Spectral Statistics Derived from Fits to Modified Periodograms. 
Periodograms Based os Prewhitened Data, and Capon s Estimators 
ofTen Synthetic 64-Point Profiles with = 1.0 • 10~ 6 and 
3 = 2.4 


Estimator 

mean c A 

std. dev. c& 

mean 0 

std. dev. 0 

Phan 

0.724 10" 6 

0.819 • 10”* 

2.767 

0.631 

Ppw 

1.310 lO" 6 

1.690 lO -6 

2.469 

0.567 

Pc AP 

0.891 lO' 6 

0.433 lO' 6 

2.371 

0.313 


Estimation of power-law parameters from the spectral 
estimates is an independent problem. In this analysis, we 
estimate c and 0 by performing a minimum absolute de- 
viation fit of a power-law function to the estimated values 
of the power density spectrum. The minimum absolute 
deviation fit, which is performed in log-log space, is more 
robust than a least-squares fit. To make the power-law fit 
independent of the frequency sampling of the spectral es- 
timates, we evaluated the estimates at frequencies spaced 
very closely together so that the estimator would be a 
smooth curve between the sampled points. While it may 
be possible to design a parameter estimate based on fewer 
samples of the spectral estimate, such a method would 
require additional assumptions in the spectral model which 
are outside the scope of this study. 

Values of the spectral estimates outside their regions of 
high accuracy (AIN < f A < 0.5 for P H AN , 1 IN < f A < 
0.2 for P pw , and 1/(2 p) < f A < 0.5 for P C ap) were dis- 
carded as before. The mean ana standard deviation of the 
estimated roughness amplitude c A and spectral slope £ are 
given in Table II. We see that the spectral parameters pre- 
dicted using Capon's estimator have roughly half the 
variance of estimates derived from either the modified pe- 
riodogram with the Hanning window or the periodogram 
based on prewhitened data. We therefore select Capon’s 
estimator as the preferred estimation algorithm for use 
with short data records in the power-law spectrum case. 

E. Formulas for Real Frequency 

Expressions for the spectral estimators listed in the pre- 
ceding sections were given as functions of the normalized 
frequency f A that is dependent on the sampling interval of 


a measured surface profile. In remote sensing studies, es- 
timates of the surface spectrum m terms of real spatial 
frequency / (in meters 1 ) are necessan so that 1) surface 
features may be compared to the electromagnetic wave- 
length and 2) spectral estimates derived from profiles with 
different sampling intervals .A ma\ be compared or com- 
bined into a composite spectrum. 

The following expressions correspond to (4). (9). ( 10). 
(6), < 12 ). ( 13), < 14). ( 15). ( 16). (2 1 ). (22), (25), and (26). 
respectively, for the case of real I nonnormalized) fre- 
quency/. where / satisfies — 1 (2A) < / < 1 (2A): 


P P fr</) = 7 ~ Z[.v„] exp < ~ / 2 7t/ A n ) 

.\ I r? - 0 

[ i:a 

EffW/H = ^ *></- $>Sz<$) ^ 

J -1 24 


"V(/> = r, 


1 / sin 7 t/A/V\ 


N \ sin 7 t/A 
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Fig. 9. Aliasing in the power-law spectrum case. Points marked (x) are 
discarded to reduce errors due to aliasing. 

F. Aliasing 

The estimates presented in the previous sections as- 
sumed that aliasing was not present, and the synthetic sur- 
face profiles were constructed to be bandlimited; i.e., they 
had no spatial frequency components above the Nyquist 
frequency f c = 1/(2A). Real surfaces will generally have 
frequency components at spatial frequencies above f c . As 
a result of a necessarily insufficient sampling rate, power 
at frequencies higher than f c will be aliased into the range 
—f c < / < f c . We can derive an equation showing this 
result in terms of real frequency by extending the deri- 
vation of Oppenheim and Schafer [23, pp. 26-28], ob- 
taining 

s. (/) = B S z (f + 2 V (40) 

In the power-law case, the power spectral density falls 
off quickly with increasing frequency. The n = 0 (exact) 
and n — 1 (first alias) terms and their sum are plotted for 
a sample power-law case in Fig. 9. We see that the effects 
of aliasing are most significant at the high-frequency end 
of the spectral estimate. As a first-order correction to avoid 
the aliasing problem, we discard values of the spectral 
estimate computed for frequencies greater than fj 2, as 
indicated in the figure. 

G. Applicability to Real Data 

Real geological data possess a variety of irregularities 
and surface statistics. Spectra of natural topography 
clearly cannot conform to a power-law model at all spatial 
frequencies— roughness features on the scale of millions 
of meters are unphysical, as are those at subatomic scales. 
However, in the studies cited in Section I as well as in 
our own investigations, measured spectra of topography 
are well modeled by a power-law spectrum in the form of 
(1) over some range of spatial frequencies. We have based 
our comparison of spectral estimators on an idealized 
spectrum with constant slope (except for the very-low- 
frequency roll-off introduced in (17) to avoid an unphys- 
ical singularity) in order to have an easy reference for 
comparison. Because both the Fourier estimators with 
windowing and Capon’s estimator are local-neighborhood 


estimators, our companson also illustrates the relative 
performance of these estimators for spectra that have 
variations in slope with frequency or local perturbations 
in level that violate monotonicity. 

III. Estimation of Two-Dimensional Spectra 

While the previous techniques allow the estimation of 
the spectrum of a linear profile, we would in general pre- 
fer a two-dimensional spectrum of surface topography. 
Such a spectrum might show hidden anisotropy or may 
indicate that the surface is isotropic in the statistical sense. 
In either case, the full 2-D spectrum provides valuable 
additional information. 

Two-dimensional spectral analyses are precluded when 
2-D data are not available. Spectral studies of seafloor 
morphology, for example, usually utilize profiles col- 
lected by bathymetric instruments that are towed by a ship, 
producing a seafloor profile along the ship’s path. Some 
directional insight can be obtained by obtaining profiles 
in orthogonal directions. 

On land, 2-D data are somewhat more accessible. 
Huang and Turcotte [1], [4] and England [5] use digital 
elevation models to estimate surface spectra, but their 
methods average over azimuth angle in the spectral do- 
main. Two-dimensional data at scales of centimeters to 
tens of meters can be tedious and costly to obtain, neces- 
sitating smaller dimensions of a sampling grid, i.e., fewer 
than 100 x 100 measured points per 2-D profile. 

A. Estimator Selection 

Two-dimensional Fourier-based spectral estimators 
such as the 2-D periodogram suffer from the same prob- 
lems as the corresponding 1-D estimators: leakage in the 
power-law case and high variance. A 2-D Capon’s esti- 
mator is one solution to these problems. The 1-D Capon’s 
estimator is readily extended to the two-dimensional case 
[20, sect. 15.8], at the cost of increased complexity. In- 
version of the covariance matrix may become problematic 
because the matrix is now of dimensions p 2 x p 2 rather 
than p x p. For a similar level of resolution, the 2-D 
Capon’s estimator will require a much greater volume of 
data (as will the other 2-D estimators). If such high quan- 
tities of data are not available, as is often the case for 2-D 
data sets, the variance of the Capon’s estimate will be 
correspondingly higher. 

B. Two-Dimensional Spectral Estimation of Isotropic 
Surfaces 

If a surface area is known to have or can be assumed 
to have isotropic statistics (from knowledge of the origin 
of the surface, or perhaps from measured 1-D spectra in 
several directions), then an estimate of the 2-D roughness 
spectrum can be obtained from estimates of the 1-D spec- 
trum. These 1-D spectral estimates may be obtained from 
individual rows of the 2-D surface height grid; we can 
calculate a spectral estimate for several or all rows and 
then average the estimates together. The resulting esti- 
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mate will have a variance that is much lower than that of 
a single row’ estimate in spite of the fact that the gnd rows 
are not completely independent. 

Suppose that an isotropic surface has a 2-D power spec- 
tral density S z ( / r , / v ) that has the form of a power law 
down to some low spatial frequency /J ovw and some un- 
specified finite form below that: 


/,./,) = Wr> = 


afr 

finite 


fr ^ /low 
fr ^ /low 


(41) 


where the radial spatial frequency / = (f \ + f 2 y ) ] We 
now need to relate the parameters a and 7 to the param- 
eters of the 1-D spectrum of surface heights measured 
along a straight line on the same surface. This 1-D spec- 
trum will also have the form of a power law for / > /j ow 
and can therefore be described by the parameters c and 0. 

To obtain this relation, we begin with the general 
expression for the Wiener-Khintchine relation in terms of 
the Fourier transform and write the correlation E[Z(x y 
y)Z(x + b x , y + 6 y )] as R Z (6 X9 by)- 


Rz(S JT> by) 



S Z (f"fy) 


• exp 0'27 t[/A + f y 5 y ]) df x df y . (42) 

Changing to radial coordinates (8 X = 6 r cos 5$, 6 V = 6 r sin 
5 e ) and noting that S z (f r9 f e ) is independent of f e , we ob- 
tain 

/* oo a 2 t 

R z (6r, 5») = I Sz(fr)fr \ exp (y2x/A cos n df df r 
Jo Jo 

(43) 

where f = / — leading to 

/% OD 

R z (6 r ) = 2x S z (f r )J 0 (2 x/ r 6 r ) / df t (44) 

Jo 


where J 0 is the zero-order Bessel function. This expres- 
sion differs from a similar (but incorrect) expression given 
by Voss [17, eq. (1.52)]. 

The 1-D and 2-D roughness spectra of an isotropic ran- 
dom field are related by an Abel transform [2]. This re- 
lation is more easily visualized using the surface autocor- 
relation functions (the inverse Fourier transforms of the 
1-D and 2-D surface spectra). The 1-D autocorrelation 
function is a slice of the 2-D autocorrelation function: 

R Z (6 X , 8 y )S(S y ) 

■ exp (-/2ir[/A + /A]) d6 x dS y . (45) 
Using Fourier transform properties: 

W/) = S z (f x ,f y ) * 6(f x ) = f s z (f x ,f y ) df y . 

J —00 

( 46 ) 



Let 5 Z | D ( /./) be given by (41). and let /, be restricted 
to the power-law region, i.e../, > / 0 »: 


W/D = 2a 



(fry' 2 d/:. 


(47) 


Changing variables, we have 


f 

Szio(fx') ~ 2a ^ f , , =• df, 

J /. v/; - /; 

which can be integrated analytically [27. p. 201 

r v ( 7 “ 1 
a vt r 

Szio(fi) ~ 


r (7/2) 


/• I -'■> - " 

Jr t 


(48) 


(49) 


where T(x) is the gamma function. 

Therefore, if S Z]D (f x ) = cf~* for f x > the 1-D and 
2-D spectral parameters are related by 

7 = 0+1 (50) 


r( 7 /2) 

c. 



(51) 


Equation (50) agrees with the result derived by Voss; he 
does not give an expression comparable to (51). 

As previously discussed, 1-D profiles of a fractal sur- 
face with 1 < Df < 2 have power-law spectra with ex- 
ponent 0 given by 0 = 5 - 2 D f9 so that 3 > 0 > 1 . From 
(50), the 2-D power-law exponent 7 satisfies 4 > 7 > 2, 
corresponding to 2-D surfaces with 2 < D f < 3, where 
D f = (8 - 7)/2. 

We next investigate the performance of the linear, av- 
eraged Capon’s estimator in the determination of the 2-D 
power-law parameters a and 7. A variation of the 2-D 
spectral synthesis algorithm described by [26] was used 
to generate 2048 x 2048-point synthetic topographic sur- 
faces with spectral exponents of 3.0, 3.4, and 3.6 with 
sampling intervals A of 1 cm. The exact surface spectrum 
was modeled as isotropic, with a form similar to (17): 


Sz(fr) = 


( (a/a 2 )|/ r | exp (-/?/[2a 2 ]) 

afP 
V. 0 


fr — /low 

/tow — fr — fc 

fr*fc 


(52) 

where f r is the radial spatial frequency, a = 2.384 • 10“ 1 1 , 
/ ow = 0.01 m -1 , and a and a are chosen as before. The 
cut-off frequency/. = 1/2 A is equal to the Nyquist fre- 
quency along the / and / axes. 

As noted previously, arrays of data containing over a 
million points are often unavailable. In practice, in situ 
surface height measurements may yield data sets that are 
much smaller, say, 40 x 40 points. Since such a data set 
is useful over a relatively narrow band of frequencies, 
several sparse grids may be collected at different scales, 
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distance, m 

Fig 10. 41 x 41 -point sample of a 2048 x 2048-point synthetic power- 
law surface with a = 2.384 ■ 10"“ and 7 = 3.0 The sampling interval 
is 50 cm. 



Fig. 11. Averaged Capon’s estimates from 2-D synthetic surface height 
data at three scales with power-law fit. 


TABLE III 

Power-Law Parameters Derived from Fits to Averaged Capon s 
Estimates at Three Scales 


Surface (exact) 

Estimated parameters j 

a 

T 

mean a 

std. dev. d 

mean 7 

std. dev. 7 

2.384 10-" 

3.0 

2.34 10-“ 

2.83 1CT 11 

2.971 

0.041 

2.384 10- 1 1 

3.4 

2.28- 10-“ 

OTiiF 

3.386 

0.041 

2.384 10-“ 

3.6 

2.17 - 10“ 11 

1.14 - 10 _li 

3.555 

0.061 


a spectral estimate calculated for each, and a function fit 
to the composite of the individual spectra. Scale factors 
are particularly important in the calculation of these 
multiscale spectra. 

In our investigation, we sample each 2048 x 2048- 
point synthetic surface at sampling intervals of 50, 10, 
and 4 cm, obtaining three 41 x 41 -point data grids. (A 
sample data grid is shown in Fig. 10.) At each scale, a 
1-D Capon’s estimate (with p = 12) is calculated for each 
grid row at frequencies V(2pA) < f < 1/(4A), and the 
estimates are averaged over the rows. A power-law func- 
tion is then fit to the averaged spectral estimates at the 
three scales. (A sample fit is shown in Fig. 11.) The upper 
frequency limit of 1/(4 A) was chosen to reduce the effect 
of aliasing. Estimates of 1-D power-law parameters c and 
(3 are then converted into estimates of a and y using (50) 
and (51). 

This process was performed on 10 synthetic surfaces at 
each of three spectral exponents. The mean and standard 
deviation of a and y are listed in Table III. 



Fig. 12. Spectral estimates calculated with (from bottom to top) P C ap< 
P H an* A>er (with no zero-padding), and Pp* from surface profiles of a de- 
bris flow near Johnston Ridge at Mount St. Helens. (Spectral estimates 
derived from techniques other than Capon’s estimator have been offset for 
clarity .) The dashed lines are power-law fits to the estimates; the equations 
of these fits are (53), (54), (55), and (56), respectively. (This and the other 
measured spectra from Mount St. Helens debris flows are presented in 
[28].). 

C. Application to Measured Data 

Surface roughness measurements were performed on 
several debris flows near Mount St. Helens in support of 
a study of electromagnetic scattering by volcanic terrains. 
The surface profiles, which were collected using several 
different sampling rates and profile lengths, were used to 
obtain estimates of the roughness spectra of the measured 
debris flows using the techniques described in the present 
report. The resulting spectra are reported elsewhere [28]. 
One of the spectral fits calculated from measurements of 
a primary debris flow near Johnston Ridge is shown in 
Fig. 12. This composite spectrum was derived from two 
surface profile grids, one with A = 1 m, the other with A 
= 1 cm. The equation of the power-law fit to the two 
Capon’s estimates is 

$z(f) = (3.57 • 10' 4 )|/r 2 - 31 . (53) 

Spectral estimates based on Phan. A>er (with no zero-pad- 
ding), and /*pw were calculated from the same data for 
comparison with (53) and are shown in Fig. 12. Equations 
of the power-law fits derived from these estimators are 

•^z,han(/) = (6.28 * 10 4 )|/| (54) 

&.PER (/) = (7.50 • 10-fil/r 2 - 36 (55) 

^z,pw (/) = (4.20 • 10 4 )|/| 24 °- (56) 

The spectral exponents in (53)-(56) seem acceptably sim- 
ilar, but this is principally due to the combination of spec- 
tral estimates at different scales into a single composite 
spectrum. If we fit power-law functions to spectral esti- 
mates based on the A = 1 cm data alone, we obtain the 
following: 

p C AP(/> = (i.i5 • io- 3 )i/r 273 


(57) 





Ifchh TRANSACTIONS UN OhN itNi'i AM) kr.Mi »! I stNMM- v < >L. ; 2 Nc i ■» K i > ^ 


v3* 


r han(/) = (3.95 • 

io- 3 )i/r ,<xj 

(58) 

PpBR(f ) = (2.12 • 

io- 4 )i/r 185 

(59) 

Ppw(/) = (9.92 • 

io - 4 )i/r : - 75 

(60) 

's, modified penodogram with Hanning 

win- 


dow. and prewhitened spectral estimators give spectral 
slopes between 2.73 and 3.0. but the penodogram (with 
no zero-padding) gives a significantly different slope 3 = 
1.85, illustrating the effects of spectral leakage. 


IV. Conclusions 

Spectral analysis is a tool of increasing importance in 
the characterization of topography and other natural sur- 
faces. Numerous studies indicate that such surfaces have 
power-law spectra in the form of (1) over some range of 
spatial frequencies. This spectral form introduces unique 
difficulties in the spectral estimation process. 

In the present work, we have shown how leakage causes 
some Fourier-based estimators to yield an estimated spec- 
tral exponent of 2.0 for surface profiles having spectral 
exponents between 2 and 3. This may explain the number 
of studies reporting power-law exponents of 2.0 for nat- 
ural terrains. We show how Capon's estimator may be 
used to avoid the leakage problem and measure the sur- 
face spectrum more accurately. We also show that 
Capon's estimator has reduced variance, which is useful 
when surface data records are short, as is often the case 
in remote-sensing studies. 

The two-dimensional Capon’s estimator may be em- 
ployed to compute 2-D spectra of rough surfaces. How- 
ever, two-dimensional data are often unavailable or too 
sparse. For surfaces that are known or can be assumed to 
be isotropic, we show how a linear, averaged Capon’s 
estimator can be used to estimate the 2-D power-law pa- 
rameters. 
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V Empirical Model 

There are no reliable scattering models for surfaces whose roughness spans 
the spatial frequencies of the radar’s wavenumber. In fact, a comparison between an 
empirical model and the predictions of several candidate theories is one of the tasks 
of the study. The empirical model was based upon scattering measurements of 
milled surfaces that were carefully produced to exhibit power-law roughness for 3 
values of B (2.29, 2.55, and 2.80). The process for manufacturing the scattering 
surfaces became Chapter V in R.T. Austin’s Ph.D. dissertation [69]. The process of 
measuring the radar backscatter from these surfaces became Chapter VI of the 
dissertation. His results were presented at an International Geoscience and Remote 
Sensing Seminar (IGARSS'94) in Pasadena, CA [68], They will also appear in a 
TGARS paper. Chapters V and VI of Dr. Austin’s thesis appears on the following 41 

pages. 




CHAPTER V 


DESIGN AND MANUFACTURE OF 
ARTIFICIAL POWER-LAW SURFACES 


The spectral estimates in the previous chapter show that certain volcanic debris 
flow surfaces at Mount St. Helens belong to the class of power-law surfaces discussed 
in Chapter II. The cited references and the present study show that power-law 
surfaces are common and often of interest in remote sensing. A better understanding 
of electromagnetic scattering by these surfaces may be obtained through scattering 
measurements and application of rough surface scattering models. 

In situ scattering measurements of natural power-law surfaces are difficult to per- 
form and even more difficult to interpret for purposes of model verification. Aside 
from common logistical problems (e.g., operation of radars in dusty environments, 
cost and time constraints for on-site work), there are several measurement uncertain- 
ties which interfere with model verification. Natural rough surfaces rarely conform 
to a sin gl e statistical description over extended areas due to the variety of processes 
that contribute to their formation. For land surfaces, volume scattering by subsur- 
face inhomogeneities is usually present. Both of these factors yield effects that are 
difficult to decouple from those due to surface roughhcfes. IKdHrl^ortnt, inter- 
pretation of scattering measurements is hindered by the difficulty of obta inin g an 
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accurate statistical description of surface roughness. As discussed in Chapters III 
and IV. estimation of the surface roughness spectrum is particularly problematic for 
power-law surfaces, requiring a large quantity of surface height data and specialized 
spectral estimation algorithms. 

The uncertainties just discussed may be largely avoided through the use of man- 
ufactured surface analogues. The analogues are homogeneous, conforming to a single 
specified statistical model. They are constructed from a homogeneous material to 
eliminate subsurface scattering, and their surface statistics are well known. Because 
the surface roughness is specified by the investigator, particular cases of roughness 
may be studied, or a single roughness parameter may be varied in order to study its 
effect on surface scattering. Finally, the scattering measurements may be performed 
under more controlled conditions, typically indoors, allowing greater accuracy and 
repeated measurements if necessary. 

This chapter describes the design and manufacture of artificial power-law sur- 
faces. The design history, material selection, and surface generation algorithm are 
discussed, and the milling system is described. The milling procedure is explained, in- 
cluding time and disk space requirements, and difficulties in milling complex surfaces 
are discussed. The chapter ends with a description of surface verification efforts. 

5.1 Analogue design history 

The original project plan called for the manufacture of power-law surfaces mea- 
suring one meter square. The surfaces were to be composed of a dielectric material 
with permittivity near that of soil. Scattering by the surfaces was to be measured 
using an X-band (9.75 GHz) radar and radars of other frequencies if available. 

Surface construction was planned to be similar to that used by Nance [44], who 
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presented scattering results from artificial Gaussian conducting surfaces [45] . As in 
his study, the surface fabrication was to be performed by a contractor; however, 
requirements for the present study were more sophisticated due to the more complex 
surface roughness (power-law instead of Gaussian) and the use of a homogeneous 
dielectric rather than structural foam that was to be coated with conducting paint. 

After an extensive search of machining companies in southeast Michigan, I was 
unable to find a contractor who was able to fabricate the surfaces for a manageable 
cost. One obstacle was the desired surface size: one meter square was too large for 
several companies’ equipment. A more significant impediment was insufficient data 
capacity. The design specification for the artificial surfaces called for the replication 
of surface features on the scale of one-tenth or one-twentieth of the radar wavelength. 
A high number of surface points were therefore required to represent the surface for 
fabrication (the surfaces eventually required 1318 x 1318 points). This quantity of 
data was quite beyond the capabilities of numerous machining firms, several of whom 
specified their data capacities in feet of paper tape. 

One company was found in Brighton which would probably have been able to 
fabricate the surfaces. However, the eventual cost per surface and inconvenience of 
working with an out-of-town company, as well as the likelihood of needing several 
attempts to get a usable surface, led to the decision to purchase milling equipment 
and manufacture the surfaces at the Radiation Laboratory. 

5.2 Dielectric materials 

Although a dielectric surface introduces greater complexity in construction and in 
backscatter modeling, the greater similarity to natural terrains led to the selection of 
a dielectric surface rather than a conducting surface. There were multiple conflicting 
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requirements for the dielectric material: 

• homogeneous (to avoid volume scattering) 

• high permittivity i to assure a measurable return signal) 

• lossy (to eliminate reflections from the bottom surface of the dielectric I 

• machinable using standard milling equipment 

• inexpensive 

• low density (to eliminate need for elaborate target mount or special handling) 

A variety of candidate materials were examined [55]. Machinable wax and mix- 
tures of machinable wax and graphite were rejected as too expensive and insuffi- 
ciently strong to support their weight (a one meter square slab might crack when 
being moved). A mixture of graphite, titanium dioxide, and silicone resin suggested 
by engineers at Texas Instruments was attempted but was never made to solidify. A 
mixture of plaster and conductive paint was rejected due to concerns about homo- 
geneity and durability (cracking under tension). The material ultimately selected, 
ultra-high molecular weight polyethylene (UHMW), was chosen for its good machin- 
ability, homogeneity, and reasonable cost. 

The search for dielectric materials led to a change in the dimensions of the arti- 
ficial surfaces. This change was principally due to mass and cost constraints. Aim 
x 1 m x 15 cm slab of these materials could measure from 130 to 260 kg, requiring 
a very strong (and probably highly scattering) surface mount and special support 
to prevent the material from cracking under its own weight. By changing the sur- 
face scale by a factor of 3.538, the X-band measurement could be simulated using a 
35 GHz radar with a smaller and shallower target surface. Surface strength was less 
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Figure 5.1: Configuration for measuring UHMW dielectric constant. 

critical at the new scale, and obtaining a measurable return was less a problem at 
35 GHz than at X-band. The decision to reduce the surface scale was fortuitous in 
view of problems with the milling table that will be described later. A surface size of 
18 inches (45.72 cm) square was selected based on availability of UHMW slabs. This 
size corresponds to an X-band surface measuring 63.7 in. (161.8 cm) square and is 
sufficiently large to allow two radar looks with minimal footprint overlap. 

Because UHMW polyethylene is not lossy, elimination of reflections from the 
bottom surface of the slabs was accomplished using the time-gating feature of the 
network analyzer that processes the radar signals. Several 6-inch slabs of UHMW 
were placed beneath the milled slab to move the bottom interface away from the 
milled surface. This configuration is described and illustrated in Chapter VI. 

The dielectric constant of the UHMW slabs was measured using the 35 GHz radar 
as shown in Figure 5.1. The radar was pointed vertically downward at a stack of 
UHMW slabs (two 6-inch and one 2-inch) and the radar bandwidth set to 3 GHz 
for maximum resolution. A metal plate was inserted at various levels in the stack to 
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identify the reflections from the top, bottom, and internal interfaces. Removing the 
plate, the time interval between the returns from the top and bottom was found to 
be 3.397 ns. The round-trip time for travel through free space at the same distance 
(0.355 m) is 2.37 ns; the ratio of these times gives the index of refraction, n = 1.433. 
The dielectric constant is then obtained directly: 

e r = n 2 = 2.05. (5.1) 

Polyethylene has extremely low loss in tabulated values below and above 35 GHz; 
it was believed that the loss tangent would be too small to measure at 35 GHz 
using techniques available in the Radiation Laboratory. The dielectric constant was 
therefore assumed to be completely real. 

5.3 Surface Design Algorithm 

The surface height was specified over a 1318 x 1318-point grid with a point 
spacing of 1.23 mm (0.04843 in.) at the X-band scale or 0.3477 mm (0.01369 in.) at 
the 35 GHz (milling) scale. The surface profile was generated using a variation of the 
two-dimensional spectral synthesis technique described by Saupe [51] and explained 
as follows: The desired surface spectrum is the two-dimensional power-law spectrum: 

Sz(fr) = a/r, (5.2) 

where f r = (fl + f^) 1 ^ 2 is the radial spatial frequency. First, a set of 2048 x 2048 
discrete Fourier amplitudes £i(/«,/ v ) is generated such that 

^|A 1 Z,(/../.)| 1 = S I (/.,/,) (5.3) 

for all points within the circle f r < / e , where / e , the maximum spatial frequency, was 
(2A) _1 = 406.5 m -1 . Point spacing in the frequency domain was A/ = (ATA)" 1 = 
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0.39698 m 1 . In other words, 

r, (/.,/,) = (5.4) 

A dimension of 2048 is chosen because the FFT algorithm requires A’ to be an integer 
power of two. and 1024 is not large enough. The amplitudes Z x are then multiplied 
by a randomness factor and a phase factor: 

Z(f t J v ) = Z l (f t ,f y )Qe (5.5) 

where Q is a Gaussian random variable with zero mean and unity variance, and 
$ is a phase random variable uniformly distributed over [0,2*]. In this way, the 
mean value of the normalized square of the spectral amplitudes still satisfies (5.3). 
To ensure that the surfac heights are real-valued, the above procedure is used to 
generate half the surface amplitudes, and the other half are obtained by enforcing 
Hermitian symmetry: 

Z(/«,/„) = £•(-/.,-/,). (5.6) 

The randomized coefficien s Z(Jm,tyi) constitute the fast Fourier transform (FFT) 
of one r eali zation of a rough surface random process having a roughness spectrum 
given by (5.2). To obtain the two-dimensional surface profile, an inverse FFT is 
performed on the coefficient array: 

Z(x k , = Z{frm, /*n)exp(— t2x[*ro + ln)/N). (5.7) 

™ msO nmO 

The resultant synthetic surface array has 2048 x 2048 points; profiles were truncated 
to 1318 x 1318 points for surfaces in the present study. 

The surface roughness parameters, a and 7 in (5.2), were selected to be s imil a r 
to those measured at Mount St. Helens and to be physically realizable on UHMW 
slabs using the milling machine. At the time the first artificial surface was designed, 
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the estimated two-dimensional roughness spectrum based upon survey site 1 and 
profilometer scans 3a and 3b had a = 1.52 x lO - * and 7 = 3.29. (The surface 
spectral estimate was subsequently changed to a = 1.7672 x lO - * and 7 = 3.345.) 
The UHMW slabs were 2 inches thick; reserving a minimum thickness of one-half 
inch for stability left a possible total relief of 1.5 inches. In order to obtain a surface 
profile having this amount of relief after scale reduction by a factor of 3.538. it was 
necessary to reduce the roughness amplitude, a, by a factor of 7. The final roughness 
spectrum used for the first surface analogue, UHMW 1 , had the mixed form of (4.50): 


where 


1 


Sz(fr) = { 


(a/<r 2 )/ r exp(-/ r 3 /[2<r 2 ]), 

«/r~^ 


f r < 0.01 m -1 
0.01 m -1 < f r < f c » 


0 , 


fr>fc 


(5.8) 


a = 2.17143 x 10' 5 
7 = 3.29 
a = 1.6438 
a = 0.00482805 
fc = 406.5 m -1 


The parameters a and a were chosen as in Chapter IV by enforcing continuity of the 
spectrum and its first derivative. In effect, however, these parameters had no effect 
(except at frequency zero) due to the limited resolution of the spectral synthesis 
technique; the lowest non-zero frequency (equal to l/NA) was 0.39698 m -1 , which 
fell outside the Rayleigh region of the mixed spectrum. 

Spectral parameters for the second and third surface analogues were selected to 
span the range of spectral slopes (or alternatively, fractal dimensions) that natural 
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Surface 

UHMW1 

UHMW2 

UHMW3 

— : '\ 

a 

2.17143 x 10" 5 ! 

2.17143 x 10" 5 

2.17143 x 10" s 

7 

3.29 

3.8 

3.545 j 

1 Q 

1.6438 

19.8523 

5.70356 

a 

0.00482805 

0.00456435 

0.00469065 

Df 

1 

2.355 

2.1 

2.2275 


Table 5.1: Spectral parameters (a, 1, a, and <r iu (5.8)) used and fractal dimension 
Dj of the surface analogues. 

terrains reasonably assume. The roughness amplitude, a, was held constant; the 
spectral slope, 7, was set to 3.8 for UHMW2 and to an intermediate value of 3.545 
for UHMW3. The values of the four spectral parameters and the associated fractal 

dimensions are summarized in Table 5.1. 

5.4 Milling System 

The milling system was ordered in November 1991 and arrived in early 1992. The 
system has existed in two forms: the original configuration, as supplied by Techno, 
and the upgraded configuration made necessary due to reliability problems with the 

Techno components. 

5.4.1 Origi al Milling System n 

The system as originally purchased (see Figure 5.2) consisted of a Techno thrae- 
axis gantry positioning table, a Kress variable-speed electronic grinder, and a Techno 
MAC 100 programmable position controller, all purchased from Techno, New Hyde 
Park, New York; CNC Software’s Mastercam CAD/CAM software, purchased from 
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Figure 5.2: The original milling system supplied by Techno. 

CIM Solutions of Canton, Michigan; and a Peer Data Systems 486/33 MHz DOS- 
compatible computer, purchased from Peer Data Systems of Ann Arbor, Michigan. 
The components of the system are described in this section. 

The three-axis gantry positioning table holds the construction material in place 
and allows the milling tool to be moved to designated xyz coordinates. The table and 
x- and z-axis supports were constructed from extruded aluminum section. Precise 
motion in the x, y, and z directions was obtained via the rotation of ball screw 
assemblies by stepper motors. The ball screw mechanisms were covered by expanding 
vinyl bellows to protect them from dust and shavings. Limit switches at one end 
of each axis’s range of motion provided a mechanism for zeroing the table. Table 
positioning was theoretically accurate to 0.01 mm (0.0004 inch), which was the linear 
motion obtained by rotating a stepper motor one step. This precision far exceeded 
the requirements for this project. The table was ordered in a special size, 48 x 
49 inches, and with extra clearance (8 inches) for the z axis to allow milling a fairly 
thick slab measuring one meter square in a single piece. 

One problem that became apparent soon after starting work with the table was 
that the table sagged somewhat under its own weight doe to support beams that were 
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insufficient for the table size. A test load of 227 kg (500 pounds) distributed over 
a 1 m 2 area in the center caused a sag of 4.3 mm. To reduce the sag to acceptable 
levels. Techno provided an iron L-beam that was attached to table surface. The 
L-beam corrected the sag problem but reduced the usable table area. Fortunately, 
the desired surface size had been reduced by this time. 

The Kress electronic grinder rotates the milling tool, removing surface material 
as it travels along a programmed path. The grinder was attached to the z-axis and 
held an endmill (milling tool) in a collet of 3/8, 1/4, or 1 /8-inch diameter. The 
rotation speed was adjustable from 8,000-20,000 rpm. 

The MAC 100 programmable position controller accepts toolpath coordinates 
from a post- processing program running on the 486 computer and translates these 
into properly timed currents in the windings of the stepper motors on the three axes. 
The MAC 100 controller was built by Techno, using some Techno-designed boards 
and other components built by Compumotor. The controller operated as an open- 
loop system; i.e., there was no position feedback from the axes to the controller, so 
the controller had no way to verify that the stepper motors and axes had moved to 
specified positions. Open-loop control is not uncommon in stepper motor systems, 
but it is less reliable than closed-loop control, as will be discussed later. 

The CNC Software Mastercam Mill package is an extensive program for designing 
complex objects, producing machine toolpaths for the fabrication of these objects, 
and controlling the machining equipment while executing the machine toolpath in- 
structions. (However, the Techno controller must be operated by a Tbchno post- 
processor program and not from within Mastercam.) Because the surfaces hr this 
study wer e extremely large and complex, aUcfc because, they were produced by || 
external program on another computing platform, many of the faattiim of ihtbtaq- 
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tercam package were not used. The principal use of the Mastercam package was to 
read a file of surface height coordinates and write files of toolpath coordinates for 
different milling tools. (Several sizes of tools were used in different stages of the 
milling process.) Techno sells a version of Mastercam designed for 286 computers; 
the faster and more sophisticated 386 version w’as needed in this study, so it was 
purchased from CIM Solutions. 

Both Mastercam and the post-processor program were run on the Peer Data 
Systems 486/33 MHz DOS-compatible computer. The computer had 16 megabytes 
of random access memory and a 210 megabyte hard disk drive, and was augmented 
with an Ethernet board to allow file transfer with the Unix workstations on campus. 
The communications link was necessary because Mastercam could process only a 
section of the artificial surface at a time. In the process, very large disk files were 
generated; it was not possible to retain all the surface files on the hard disk. 

5.4.2 Failure of the Techno System 

The Techno system began having problems in October 1992. The z-axis motor 
either slipped or rotated in the wrong direction and seemed to be overheating during 
some test runs on wood surfaces. Upon inspection, the solder joints on the connector 
card inside the motor housing had heated up enough to melt the solder. Engineers 
at Techno were unable to determine a cause for the failure. 

In January 1993, it was proposed that the z-axis motor might be overheating due 
to excess weight on the z axis. The z axis was reverse-mounted on the Techno table 
to allow the full 8 inches of clearance; i.e., it was installed so that the movable carrier 
stayed fixed while the rest of the z-axis assembly moved up and down. To relieve 
part of the load on the z-axis motor, an 11.5-inch tension spring and support rod 
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endmill 


Figure 5.3: A load- bearing spring was installed on the z-axis to help support its 
weight. 

were installed as illustrated in Figure 5.3. The spring constant and position were 
selected such that the spring force balanced the weight of the z-axis assembly in its 
mid-position. The spring seemed to solve the slippage problems. Milling on the first 

UHMW slab commenced in February 1993. 

By early March, the machine seemed sufficiently reliable to leave it operating 
overnight, which was a great advantage, considering the very lengthy milling times. 
However, the machine began to malfunction again in late March: this time, the y 
axis became confused during a finishing run, resulting in a deep gouge in the slab. 
Manual commands to the mill were sporadically executed incorrectly. At this point, 
AC line noise was suspected, so an isolation transformer was installed. Machine 
operation was then normal again. 

In mid- April, the problem recurred. Both the y and z axes malfunctioned during 
an overnight roughing operation. The 1/4-inch endmill had milled doWn through the 
slab and through part of the table itself. The simultaneous failures of two axes led 
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Figure 5.4: The upgraded mill configuration. 

to the conclusion that the Techno MAC 100 controller was not reliable. This was 
a reasonable assumption in view of the other engineering problems with the Techno 
equipment. 


5.4.3 Upgraded Milling System 

The upgraded milling system is shown in Figure 5.4. The Techno MAC 100 was 
replaced by a Compumotor AT 6400 Indexer and a Compumotor SC30 Motherboard 
Rack, purchased from Amerinetics of Novi, Michigan. Optical encoders (part num- 
ber MOD5641-25-100-L) were purchased from BEI Motion Systems of Plymouth, 

, *i ' t 

Michigan. Motion Architect software, included with the AT6400 from Compumotor, 

. . : • .0- - <s 

was used with Microsoft Visual Basic to write new control software. The indexer. 
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drive rack, and encoders arrived in late July 1993. 

The AT6400 indexer is the central controller of the upgraded system. The indexer 
receives ASCII-format move instructions from the control software, for example, 
“move x axis 100 steps right.” The indexer then converts these instructions to step 
pulses and direction signals that are sent to the stepper motor drive card for each 
axis. These drive cards were removed from the MAC 100 but were manufactured by 
Compumotor; they were installed in the Compumotor SC30 Motherboard Rack for 
use in the upgraded system. The three stepper motor drivers use the step pulses and 
direction signals to drive the windings in the stepper motors in sequence, producing 
rotation in the specified direction. 

The optical encoders were installed to provide feedback to the indexer. They 
were mounted on the stepper motor cases using epoxy such that the motor shafts 
projected through the encoders. The encoders consist of two LED’s and a rotating 
glass disc with 100 equally spaced lines. As the encoder turns, light from an LED is 
alternately blocked and then transmitted through the glass disc. The two detectors 
are in quadrature, resulting in 400 encoder pulses per revolution. The encoders are 
relative sensors — they do not provide absolute coordinates, so the indexer must keep 
a count of encoder pulses. The control software was written so that after each mill 
move, the encoder step count was compared to the stepper motor step count. If the 
counts disagreed by more than two steps, the program stopped, and the mill was 
deactivated. 

In addition to the encoder inputs, the indexer also had inputs from a home limit 
switch on each axis and a limit switch on the z axis that stopped the mill if the z 
axis went too low. Additional inputs from jog switches allowed manual control of the 
mill position. The indexer also had programmable outputs, one of which was used 
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to operate a relay that activated and deactivated the grinder motor. 

One other component of the MAC 100 was retained: the large transformer that 
converted 115 VAC down to 18 and 26 VAC to supply the motherboard rack and 
stepper motor drives. A transformer malfunction seemed unlikely. 

5.5 Milling procedure 

5.5.1 Generation of Toolpath Files 

Surface fabrication begins with the spectral synthesis procedure, described earlier 
in this chapter, which results in an 22.6 MB ASCII file of surface elevations arranged 
in a 1318 x 1318-point grid. Because the Master cam software cannot accept input 
of more than 100 lines of ASCII surface elevations, the surface must be segmented 
into thirteen sections of 99 lines plus a fourteenth section of 31 fines. The FOR- 
TRAN program ds2mcinv .f performs this segmentation and applies the scale factor 
converting the surface from the X-band scale down to the 35 GHz scale. The pro- 
gram writes the surface sections as long fists of xyz coordinates. The resultant files 

measure about 5 MB each. 

The edges of the surface analogues were rounded by subtracting the random 
scaled height from a 1 /4-inch radius along each edge, as shown in Figure 5.5. The 
rounding was performed so that no non-random sharp edges would be illuminated 
by the radar. Most rays striking the rounded edges should be scattered away from 
the radar. The FORTRAN program adganax.f reads the surface file and writes a 
3.8 MB file of edge coordinates. 

A couple of sections at a time are downloaded to the 486 computer. Mastercam is 
then used to convert the surface elevation files into toolpath files at two mentions. 
Roughing files contain toolpaths for a 1/4-inch ball-end endmffl. This tool is used to 
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Figure 5.5: Edges of the artificial surface are rounded by subtracting the surface 
elevation from a rounded reference surface with a 1/4-inch radius of cur- 
vature. 

remove the bulk of material at a rate much faster than the smaller tools. Finishing 
toolpaths are much more detailed (and much larger) and are used for the passes 
using the 1/16-inch and 1 /32-inch ball-end endmills. Cross toolpaths are for the 
orthogonal finishing pass with the 1/32- inch endmill. 

The toolpath ( .NCI) files are produced in several steps [54]. An ASCII file is input, 
and the resultant surface description is saved in Mastercam format as a geometry 
(.GE3) file. The next step is to chain the surface, i.e., to define the contours to 
be milled. This creates a chaining (.NCS) file. After chaining, which need be done 
only once per section, the tool parameters and desired resolution are entered, and 
Mastercam writes surface (.CDB), offset surface (P.CDB), and toolpath (.NCI) files, 
in addition to some miscellaneous small intermediate files. When all three toolpaths 
have been written, the roughing toolpath is filtered by removing points within a 
specified distance from the splines used to define the surface. This significantly 
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File type 

extension 

roughing 

finishing 

cross 

geometry 

* . GE3 

12 

chaining 

*.NCS 

12.5 

12.5 

12.5 

surface 

*.CDB 

0.4 

3.2 

3.2 

offset 

•P.CDB 

0.4 

3.2 

3.2 

toolpath 

* .NCI 

0.02 

4.6 

4.6 


Table 5.2: Types and sizes (in MB) of Mastercam files. 

decreases the size and execution time of the roughing toolp&th. 

Entering, processing, and transferring the generated files for a single section takes 
about 3 hours. The sizes of the various Mastercam files are given in Table 5.2. Some 
of these may be deleted (the offset surface and all but one of the chaining files), but 
the rest are uploaded to a Unix workstation and saved on magnetic tape, to avoid 
lengthy re- processing if a toolpath file is lost. 

5.5.2 Milling machine operation 

The surface slab, which has been milled fiat using a flat-end endmill, is mounted 
on the milling table after aligning its edges with the x and y table directions. It is held 
in place using aluminum angles tightened against each side of the slab. The mounting 
brackets are tightened carefully because re-aligning the slab would be extremely 
difficult if it came loose. In several instances, it was necessary to mill the ahiminnm 
brackets down a few millimeters in places where the surface was low near its edge. 

A Visual Basic program called mill controlled the entire milling process. Vi- 
sual Basic was quite useful because it allowed easy construction of virtual control 
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panels and readouts describing the milling machine status. During operation, the 
program displays the current position of the mill (in steps and in inches), the encoder- 
determined position (in steps and in inches), and an estimate of time remaining until 
completion. 

To begin roughing, a 1 /4-inch ball-end endmill is placed in the grinder collet and 
tightened firmly. After homing the mill to zero using the limit switches on the three 
axes, the endmill is positioned over the starting corner of the surface. It is lowered 
carefully until a sheet of paper can no longer slide between the endmill and the 
surface. This position is recorded as the reference level. Next, the endmill is moved 
over a small block of metal attached to the table near the slab, and the process 
is repeated. The z coordinate is recorded and is used later to verify the endmill 
position. For example, the indexer and program have no way to detect when an 
endmill slips within the collet, but the metal block test will show the slippage. The 
block also provides a way to determine the starting position after changing tools. 

The first roughing pass is offset above the reference plane so that the maximum 
depth milled is 0.2 inches. Subsequent passes remove 0.2-inch layers; roughing is 
continued until only 0.1 inch of material to be removed remains. The electronic 
grinder is set to speed 2 for roughing. Roughing takes about 45 minutes per section 
per pass. 

An intermediate pass is next made using the finishing toolpath and a 1/ 16-inch 
ball-end endmill. The intermediate pass is offset 0.025 inches above the final surface. 
The grinder is set to speed 1 for the intermediate and finishing passes. The slower 
speed is necessary due to the reduced rate at which shavings are cleared from the 
flutes in the smaller endmills. (Shaving buildup causes the endmill to become hot.) 
The intermediate passes last about 5 hours per section. 
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radius = 0.0156 inch 



Figure 5.6: The spacing between adjacent finishing rows (0.0137 inch) was slightly 
less than the radius of the 1/32-inch ball-end endmill (0.0156 inches). 

The edges are milled after the intermediate pass using the 1/ 16-inch endmill and 
the toolpath written for the edges. Two edge passes are completed: one with offset 
0.1 inch, and another with an offset of 0.025 inches. 

The finishing passes are run using a 1/32-inch ball-end endmill. The miniature 
endmills were purchased from R & F Micro Tool Co. of Pembroke, Massachusetts. 
Finishing passes last 5 hours as well. In milling test wood surfaces and the first 
UHMW surface, it was determined that often the surface material was not completely 
removed in the finishing passes. This was partially due to the limited overlap between 
adjacent rows when using the 1/32-inch ball-end tools (as shown in Figure 5.6), and 
partially due to the endmills becoming less sharp with use. The residue remaining 
could sometimes be pulled off in strips, but often it formed individual fibers that could 
be described as plastic “peach fuzz.” Sanding, singeing, and brushing were tried in 
attempts to remove the fibers, but the procedure found to produce the fewest fibers 
was to make the initial 1/32-inch finishing pass in the cross direction (i.e., orthogonal 
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to the section’s long dimension) using zero offset, followed by a lengthwise finishing 
pass (i.e., along the section) with a positive offset of 0.002 inches. This procedure 
reduced the fibers to acceptable levels. 

Keeping the tools sharp was a constant task. The endmills tended to become dull 
quickly due either to the material or to the high rotation speed. Tool dullness was 
particularly troublesome with the 1 /32-inch endmills, which have very little cutting 
area, and which need to be the most sharp. Some of these endmills seemed dull when 
new. After consulting with R k F Micro Tool, endmills coated with titanium nitride 
(TiN) were found to keep a sharp edge longer. 

The finishing process was quite lengthy, with 28 section passes at 5 hours each. 
The mill ran unattended, but was checked by an operator every few hours. Three 
sections per day was the effective maximum rate of completion. 

5.6 Surface Completion and Temperature Sensitivity 

Programming and testing of the upgraded system continued through September 
1993. Im provements during this period included the addition of a blower fan to keep 
the z-axis motor cool and optimization of mill travel speeds for roughing, intermedi- 
ate, and finishing passes. In late September, work was started on the first artificial 
surface, UHMW1. 

A new problem appeared during the intermediate-scale milling on UHMWl. Af- 
ter noticing that the slab was slightly loose within the holding clamps early in the 
morning, a heat gun was used to heat the surface slightly. The slab soon became 
immobile, having expanded due to the increase in temperature. The thermal expan- 
sion and contraction became a significant problem in the milling process. When the 
temperature dropped, the slab became loose; when the temperature rose too high, 
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the slab bowed upward because the mounting brackets prevented lateral expansion. 
Either of these changes left visible ridges in the milled surface. The thermal contrac- 
tion and expansion tended to lag behind the room temperature by several hours. A 
heat lamp and space heater were tried as means of keeping the surface warm, but 
neither could keep the slab at a constant temperature. In October, the steam supply 
to the room heater was turned on, but the radiator was unable to keep the room 
above 20 °C at night even at its maximum setting. (The milling machine was located 
in room 422-4 in the Aerospace Engineering Building; this room is quite large with 
a 20-foot ceiling and has windows on three sides which are poorly sealed.) 

The ambient temperature had been 24-27 °C when the surface was first mounted; 
maintaining this temperature became more difficult as outside temperatures fell 
through October. After repeated maintenance work on the room heater and im- 
proved sealing on the windows failed to help, a cubical enclosure (“the Greenhouse”) 
was constructed from light wood and plastic sheeting around the milling machine. 
A compact electric heater with thermostat and a blower fan were placed inside in 
an attempt to stabilize the the temperature. The electric heater could not maintain 
an 80 °F air temperature around the clock, but it did provide sufficient heat during 
the day to allow UHMW1 to be completed on 19 November. A verification scan was 
performed on the surface before removing it from the mounting clamps. The scan is 
described in the next section. 

Milling on UHMW2 began on 11 March 1994 and was completed on 14 April. 
Temperature stability was less problematic because the milling was started at a lower 
temperature that was maintainable by the electric heater. 

The milling machine was moved to room 422-15 before starting UHMW3 on 
14 June. Room 422-15 had an air conditioner, which was necessary to avoid slab 
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expansion due to overheating. UHMW3 was completed on 9 July. 

5.7 Surface Verification 

After each surface was completed, the electronic grinder was removed and a small 
ranging laser mounted on the milling machine. The new laser, which was purchased 
to replace the surveying laser used in Chapter III, produced elevation measurements 
over a limited range at a very rapid rate. The laser generated an output voltage 
proportional to the measured distance; this voltage was read by a sampling multime- 
ter and relayed to the 486 computer through an HP-IB link. The surface scanning 
program, nillscan, was written in Visual Basic and was quite similar to the milling 
program. 

The laser beam was carefully adjusted to be vertical and was aligned as closely as 
possible to the origin of the milled surface. The surface scan consisted of 264 x 264 
points with a point spacing equal to five times the point spacing used in the original 
surface specification. The elevation at each point was based on the average of 100 
samples measured at two-millisecond intervals. 

An error estimate was calculated by subtracting the scanned elevations from the 
elevations specified in the original surface file. The error estimate from UHMW1 is 
shown in grayscale in Figure 5.7. No systematic features are visible in the file. The 
surface in the figure is actually a subset of the error estimate; several rows around 
the edge were removed because the laser range estimates were visibly wrong near 
the mounting brackets. The corrupted readings were most likely due to the angled 
orientation of the detector which senses the reflected laser light. This configuration 
seemed to decrease the laser accuracy on sloped surfaces, suggesting that the width of 
the error distribution (shown in Figure 5.8) was inaccurately large. The repeatability 
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Figure 5.7: Grayscale file of estimated error in UHMW1, calculated by subtracting 
the scanned height from the specified height in the original surface file. 



surface error (inches) 

Figure 5.8: A histogram of differences between the specified surface elevations and 
the scanned surface height. 
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of the endmill position on the metal block position to within 3-6 steps (0.03-0.06 mm) 
after numerous milling sessions supports the hypothesis that errors in the milled 
surface elevations are below the resolution of the scanning laser. 

Fabrication of the surface analogues was by far the most time-consuming part of 
this study. Measurements of radar backscatter from these surfaces is described in 
the next chapter. 




CHAPTER VI 


RADAR MEASUREMENTS 


Backscatter from the artificial power-law surfaces was measured at 35 GHz us- 
ing the Radiation Laboratory Millimeter-wave Polarimeter. Measured values of the 
scattering coefficients of the scaled surfaces were then interpreted as X-band scatter- 
ing coefficients of equivalent larger-scale surfaces. (The scale factor turns out to be 
unity.) 

In this chapter, the radar components and measurement configuration are de- 
scribed. The calibration and measurement algorithms are discussed, and final post- 
processing and data reduction steps sue detailed. 

6.1 35 GHz Radar 

The Radiation Laboratory Millimeter-wave Polarimeter (MMP) was first assem- 
bled in 1986 [62] and was most recently rebuilt in 1994. The MMP is a dual-antenna, 
network- analyzer- based scatterometer designed for operation at 35 GHz. A phase- 
locked on-board local oscillator signal is mixed with the 2-4 GHz IF sweep received 
through cables from the network analyzer to produce a transmitted wave sweep from 
34 to 36 GHz. The received signal is downconverted by mixing with the signal from 
a second local oscillator, resulting in a 2-4 GHz IF signal that is transmitted to the 
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network analyzer through coaxial cables. The transmit and receive local oscillators 
are phase-locked by split signals produced by am on-board X-band oscillator that are 
multiplied up to the LO frequency. Motorized rotating wave plates in the transmit 
channel allow the transmission of arbitrary polarizations, and separate horizontal (H) 
and vertical (V) receive channels allow the simultaneous measurement of H and V 
return signals when using a network analyzer with dual-channel reception capability. 
The transmit and receive assemblies may be separated for use as a bistatic radar. In 
addition, the removable transmitter module may be installed in the receiver assembly 
allowing monostatic operation using a single antenna. A block diagram of the radar 
is shown in Figure 6.1. 

For this study, the radar was configured as a dual-antenna, quasi-monostatic 
instrument operating at the near edge of the far-field region. The transmit and 
receive antennas were installed adjacent to each other and offset by 4.8° so that the 
transmit and receive antenna boresights converge at the target range. The 3 dB 
product beamwidth was 3° and the 6 dB product beam width measured 4.2°. While 
the original radar specifications called for components with a bandwidth of 2 GHz, a 
slightly wider bandwidth of 2.4 GHz (33.6-36 GHz) was used for the measurements 
in this study. 

An HP 8720 network analyzer was used to process the returned signals and trans- 
late them for output to a Gateway 2000 DOS-compatible computer. The HP 8720 
was preferred over an HP 8753 because its frequency range allows input and output 
of signals up to 4 GHz, making IF sweeps of 2-4 GHz (or 1.6-4 GHz) possible. This 
advantage was offset by the test set of the HP 8720, which has only two ports, pre- 
cluding reception of both V and H channels simultaneously. An external microwave 
switch was used to connect either the V or H receive channel to the network analyzer. 
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This switch was operated by the computer using an HP 3488 A switch/control unit. 

The time-gating capability of the network analyzer allows the operator to isolate 
the response from a selected range interval in the time domain. This feature was used 
in the present study to mask out reflected signals from the bottom of the UHMW 
stack, the turntable, the wall behind the target, and other objects in the room, 
making accurate cross section measurements possible without an anechoic chamber. 

The network analyzer, transmit polarization switches, receive channel switch, and 
target turntable were all controlled by a Gateway 2000 4DX2-66V DOS-compatible 
computer using control software written in Visual Basic. Raw and time-gated trace 
data were stored on the computer’s hard disk drive. The radar and measurement 
control system is shown in Figure 6.2. 

6.2 Measurement Configuration 

The scattering measurements were performed in a side area of the Radiation Lab- 
oratory Bistatic Scattering Facility, Room 422-4, Aerospace Engineering Building. 
Figure 6.3 shows a measurement schematic. 

The milled surface was placed atop an 18-inch stack of UHMW slabs (two 6- 
inch slabs on top of three 2-inch slabs) in order to increase the delay between the 
reflections from the top and bottom surfaces, allowing use of a time gate to isolate 
the upper surface reflection. The slabs were flat except for a small air gap around 
the edges of the two 6-inch slabs and a slightly larger air gap (maximum thickness: 
1 mm) between the milled surface slab and the top of the UHMW stack. This 
larger air gap was due to a slight warping of the milled slab caused by the release 
of internal stresses during the milling process. A dielectric filler (Vaseline petroleum 
jelly, tr = 2.16 — ./ 0.0022 at 10 GHz [29]) was used to fill the air gaps. Plots of 
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Figure 6.2: Components of the radar measurement system. 
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Figure 6.3: Schematic of measurement setup. 

the time-domain response before and after applying the filler showed a reduction in 
reflections from the internal interfaces. 

The UHMW stack was placed on a square piece of Eccosorb AN-PXP-74 foam 
absorber. Surfaces UHMW1 and UHMW2 were measured with the slab stack and 
absorber placed on the large turntable used for the sandbox in the bistatic system. A 
much smaller turntable was used for the UHMW3 measurement; a square of 1 /4-inch 
plexiglass was placed between the turntable and absorber to increase stack stability. 

The target turntable was placed about a meter in front of an absorber wall to 
eliminate reflections from the corner of the room. Other objects were removed so 
that the artificial surface was the sole object at target range near the main beam of 
the antenna. A large section of absorber was placed in possible sidelobe directions, 
but no significant reflections due to sidelobes were detectable. 

The radar transmitter and receiver assemblies rotated about a front axis that 
was mounted on the front of a large steel frame. The incidence angle was adjusted 
by pivoting the radar around the front axis, using a small electric motor to reel in a 
cable attached to the rear plate of the receiver. The steel frame was attached to a 
manually operated, mechanically telescoping lift (a Genie Superlift, manufactured by 
Genie Industries of Redmond, Washington). The lift was used to raise the radar to 





Ill 


heights up to 3.5 m, making measurements at near-normal incidence angles possible. 

The system had a high signal-to-noise ratio. The signal peak corresponding to a 
1.75-inch sphere was about 58 dB above the noise level before calibration. Isolation 
between the co- and cross-polarized channels was about 25 dB before calibration and 
about 42 dB after calibration over most of the bandwidth. 

6.3 Calibration 

The radar was calibrated using a recent technique for polarimetric coherent-on- 
receive radar systems [46]. The technique uses two calibration targets (a metallic 
sphere and a depolarizing target) for the initial system calibration and a sphere 
alone for subsequent gain corrections. A 1.75-inch (4.445 cm) diameter steel sphere 
and a pair of 0.028 cm diameter x 5.4 cm lengths of wire mounted in styrofoam were 
used as the calibration targets for these measurements. The calibration technique is 
described in Appendix B. 

The radar system stability needed to be established before measurements of scat- 
tering by the artificial surfaces could begin. Repeated calibrations using a sphere and 
depolarizing target were used to test system stability as reflected in the calibration 
coefficients. Several problems were fixed in order to get repeatable values for the 
time-invariant calibration constants Ci, ca, C 3 , r lt and ra: 

• The use of trace subtraction to remove the response from the calibration target 
mount was discontinued because the mount position changed when placing or 
removing the calibration target. Fortunately, the scattering cross sections of 
the calibration targets were much higher than that of the styrofoam mount. 

• A sturdier styrofoam column was used. The steel sphere used as a calibration 
target was heavy enough to compress and give a slight tilt to the original 




styrofoam mount. 


• A phase-lock problem was fixed by slightly changing the frequency of the X- 
band oscillator used to phase lock the two local oscillators. 

• The radar temperature was stabilized, first by installing sheets of clear plastic 
to make the radar assembly near-airtight, and later by installing a heater inside 
the enclosure. 

• The calibration target was allowed to settle for at least 15 minutes after placing 
it on the styrofoam target mount. Repeated tests indicated a non-visible sway 
introduced by placing the target (especially the heavy sphere) on even the 
sturdier styrofoam column. Target stability is particularly important while 
measuring ci, cj, C 3 , T \ , and tj; the settling period increased the repeatability 
of these factors. 

• A minimum warm-up period of one hour was selected for both the radar and 
network analyzer. 

These measures, together with more mundane practices like avoiding contact with 
the radar cables or mount, brought the system to a sufficient level of stability, as 
reflected in the time- invariant constants and the channel gains R\ and R-j. 

Test measurements of the calibration sphere (Figure 6.4) and a smaller sphere 
measuring 1.188 inches in diameter(Figure 6.5) showed that deviations from theoret- 
ical scattering cross sections were less than 0.5 dB (usually less than 0.3 dB) and that 
cross-polarized isolation was good, usually 40-45 dB or more. Although the system 
stability seemed good over fairly short periods, it was suspected that longer-term 
variations would be greater. A conservative approach of calibrating before and after 
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Figure 6.5: Test measurement of 1.1 86- inch diameter sphere. 
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each rough surface data set was adopted. During the processing of the measured 
data from UHMW1 and UHMW2, it was noticed that the magnitudes of /?, and 
R 2 computed after a measurement run were slightly different from those measured 
before the run. The relative phases also varied as a function of frequency. Linear 
interpolation of the magnitudes and phases of Ri and R 2 with time was used as a 
first-order correction to minimize the effect of the changing R parameters. Values 
of Ri and R 2 for a given trace are based on linear interpolations computed with 
respect to the actual times of the previous and subsequent calibrations and the start 
and stop times of the measurement set. Test measurements in which the phases 
of Ri and R 2 were 180° off showed that errors in a were only about 0.6 dB; the 
linear-interpolation-based values should be considerably more accurate in phase. 

6.4 Measurement Procedure 

After characterizing the system by performing a full calibration to determine 
the time-invariant calibration coefficients, the UHMW stack was placed upon the 
turntable, and the surface to be measured was placed on the stack, as discussed in 
Section 6.2. Dielectric gel was used to fill in air gaps under the edges of the surface 
to be measured. 

Before beginning measurement, the calibration sphere (1.75 inches in diameter) 
was placed on its mount. The radar was moved to the operating range for this 
study, 2.7 m. (The range was measured using a 2.7-meter length of wire attached 
to the radar faceplate.) After carefully peaking the sphere signal by adjusting the 
radar azimuth and elevation angles, a red targeting laser (attached below the receiver 
antenna) was adjusted to point at the center of the antenna beam at the operating 
range. The laser made later pointing much simpler. 
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The artificial surfaces were measured at incidence angles ranging from 15° to 60° 
at 5° increments. The measurement sequence at each incidence angle consisted of 
the following steps: 

• with radar lowered, set incidence angle, 6 

• set h, radar height, and d, distance from lift base to target surface 

• raise radar in elevation to point at sphere and calibrate 

• reset incidence angle, 0 

• raise back to height h and position using target laser and range wires 

• measure 60 samples of surface backscatter using W, HV, VH, and HH polar- 
izations, rotating the target 5° between measurements 

• raise radar in elevation to point at sphere and recalibrate 

A few additional sphere measurements were inserted between calibration sets. 

6.5 Post-Processing and Data Reduction 

During the measurement runs, raw trace data is saved to disk files in binary form 
for reasons of speed (it is much faster than saving gated data) and flexibility (the 
data may be gated later in various ways). 

Gating refers to the process of transforming a frequency-domain data trace to the 
time domain; applying a time gate (i.e., a bandpass filter of specified width that is 
used to eliminate signals from objects at ranges outside the gate), and transforming 
the result back to the frequency domain. The network analyser was set to its standard 
gate shape; the gate width was set to 5 ns in order to indude reflections frocp both 
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near and far parts of the illuminated surface while excluding reflections from the 
bottom of the stack. The gate width was reduced to 4 ns for incidence angles of 25°, 
20°, and 15° due to the smaller range difference between the closest and most distant 
surface points. The gate center was set at the mean target range, usually 18.01 ns. 

After the surface and calibration target files have been gated, they are uploaded 
to a Unix workstation for further processing. Measurements of the calibration target 
are first converted into values of R\ and flj as a function of frequency (i.e., for each of 
the 401 discrete frequencies that make up a data trace). Next, files of gated surface 
data are converted into average radar cross sections, <x, for each polarization as if the 
targets were point targets (rather than distributed targets) using 

»,,(/.) = <IV^«(/,)| J ), (6-D 

where r and t indicate the receive and transmit polarizations and may be either V 
or H, and 5 rt is an element of the scattering matrix given by (B.43). The angle 
brackets {•) in (6.1) denote averaging over all traces in a set (i.e., over all surface 
measurements at various azimuth angles). Averaging is done separately for each of 
the 401 trace frequencies. 

6.5.1 Scattering Coefficients 

Values of the scattering coefficient er°(/,) are obtained by applying an illumination 
integral correction to cr(/,). The correction may be derived from the forms of the 
radar equation used for point and distributed targets. The point target form for a 
target at boresight is [62] 

*(«) = f,< ^T A * jF g W- (6 - 2) 




117 


tu 


and the distributed target form is 


Pr(6) = 


PtGot Gor 
(4r) 3 


I(h.B)o°(9). 


(6.3) 


where I(h.6) is called the illumination integral: 


'< M > = /„ 


gt(9ti 0t)9r{^r> 4>r) r^, n 
Ilium, axes 0) 


dA. 


(6.4) 


The illumination integral accounts for the variation in antenna gain and range 
across the distributed target. The illumination integral correction is obtained by 
setting the ratio P r {0)/P t equal in (6.2) and (6.3) and solving for <7° in terms of <r: 

'°<*> = (6 ' 5) 

The program measillum.f was used to convert values of <r(/,) to values of <7°(/,) 
for four polarizations at 401 frequencies per data set based upon the range and 
incidence angle for each measurement. Normalized antenna patterns were based 
upon data supplied by the manufacturer. 

Final values of < 7 ° were obtained by averaging over frequency. Errors in the 
measured cross section of test spheres were greatest at the upper and lower ends 
of the frequency band; frequencies outside the range 33.8-35.8 GHz were therefore 
omitted from the frequency averaging. 

The values of the scattering coefficient were obtained using a surface model scaled 
for use with a 35 GHz radar instead of an X-band radar. Using Table 11-6 of 
Ruck [50], we see that the scaled quantities (denoted by primes) of length, conduc- 
tivity, frequency, wavelength, permittivity, permeability, and radar cross section are 
related to the unsealed quantities as follows: 


/' = IfP, 


( 6 . 6 ) 
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pa, (conductivity) 

(6.7) 

p/. 

(6.8) 

Vp> 

(6.9) 


(6.10) 

Pi 

(6.11) 

er/p 2 , (radar cross section) 

(6.12) 


where p is the scale factor (3.538 in this study). The conductivity scaling is not a 
problem because UHMW has negligible conductivity. 

The scale relation for the scattering coefficient of a distributed target, <r°, is 
obtained by noting that the surface area, A, scales as length squared: 

*=<o’ = g)’ = £. < 613 > 

Writing the scattering coefficient as the average radar cross section normalized by 
area, we obtain 

(„0 y _ W = Wifi = M = „ 0. (6.14) 

1 ' A' A/r> A 

Because <r° is a normalized quantity, it is independent of the model scale. 


6.5.2 Independent Samples 

Independent samples of the surface scattering were obtained through a combi- 
nation of target rotation, spatial sampling, and frequency averaging. The radar 
6-dB footprint was positioned as close to the surface edge as possible to minimize 
overlap between footprints (Figure 6.6). Overlap of the 3-dB footprints was smaller 
still. The surface was rotated 5° between measurements. A few measurement runs 
were executed using smaller angular increments, A <f>, in order to study the decorre- 
lation as a function of target rotation. Based on limited angular decorrelation sets, 




6 dB 

footprint 

Figure 6.6: Arrangement of footprints on slab 

decorrelation angles for UHMW3 are approximately 1.1° for 60° incidence and 4.5° 
for 15° incidence (based on 60-sample traces with A^ = 0.5°). For UHMW1, the 
decorrelation angle is approximately 1.5° for 60° incidence (based on a 90-sample 
trace with A <f> = 2°). A data set with fine angular increments was not taken on 
UHMW1 at 15° incidence, but estimates based on the regular data sets indicate a 
higher decorrelation angle, though probably still less than 10* . No decorrelation sets 
were collected on UHMW2. The number of independent samples obtained through 
surface sampling N„ therefore varied from 300*/7.5* = 40 (estimate for UHMW1 at 
15°) to 60 (the number of samples collected) for UHMWl at 60° and UHMW3 at all 
incidence angles. 

Values of < 7 ° were averaged over a 2 GHz bandwidth, increasing the effective 
number of independent samples by a factor dependent on the footprint’s extent in 
range. Ulaby et al. [57] give an expression for the effective number of independent 
samples obtained through continuous integration over a swept-frequency bandwidth 
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surface 
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N„ 

— 

Njavg 

Ntotai 

UHMW1 

15“ 

40 (est.) 

1.3 

52 


O 

O 

ID 

60 

5.1 

306 

UHMW3 

15° 

60 

1.3 

78 


60° 

60 

5.1 

306 


Table 6.1: Estimated number of independent samples 


B: 



(6.15) 


where a = 2xD/c and D is the difference in range between the nearest and most 
distant points in the radar footprint. Using the dimensions and bandwidth from the 
present study, Nj avg can be shown to vary from 1.3 at 15° incidence to 5.1 at 60° 
incidence. 

The total number of independent samples, N to tai, is calculated by multiplying the 
numbers gained by spatial sampling and frequency averaging: 


Ntotai N tt Nf a vg 


(6.16) 


The resulting numbers are listed in Table 6.1. Values of N to tai for UHMW2 were 
probably slightly lower than those of UHMW3 due to decreased roughness. 

Values of the backscattering coefficient measured for the three surfaces are given 
in Chapter VII. 
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VI A Comparison with Various Scattering Models 
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under' this investigation, the Smail Perturbation 

Model performed best. 

The comparison among theories and experimwual 

on the following 14 


pages. 




CHAPTER VII 


ROUGH SURFACE SCATTERING MODELS 
AND COMPARISONS TO EXPERIMENT 


The surface scattering measurements described in the previous chapter were per- 
formed in May and July of 1994. In this chapter, results of those measurements are 
shown and compared to predicted values of a° resulting from several rough surface 
scattering models. 

7.1 Experimental Results 

Measured values of and <tJJa for the three artificial surfaces are shown in Fig- 
ure 7.1. The open symbols represent W polarization and the filled symbols represent 
HH. Backscatter dependence on the surface roughness behaves as expected, with the 
roughest surface (UHMW1) having the highest backscatter, followed by the interme- 
diate surface (UHMW3) and the smoothest surface (UHMW2). The backscattering 
coefficient decreases with increasing incidence angle for all three surfaces, as ex- 
pected, although there is a small upturn in a° values for both polarizations at 5.5° 
for the two smoother surfaces. There is little difference between the and 
for a given surface — is usually slightly higher, but never more than a couple of 
decibels. The difference becomes small at small incidence angles, as it should. 
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incidence angle (degrees) 

Figure 7.1: Measured co polarized backscattering coefficients, and <r° A , of the 
three artificial surfaces: Surface A (UHMWl), Dj = 2.355; Surface B 
(UHMW3), Dj = 2.2275; and Surface C (UHMW2), Dj = 2.1. 

Measured values of the cross-polarized backscatter, <r° h and <r® u , are shown in Fig- 
ure 7.2. The measured <r° v and also follow expected trends, decreasing with in- 
creasing incidence angle, and decreasing for smoother surfaces. The measured values 
are generally 20-23 dB lower than the co- polarized backscattering coefficients. Values 
of and <7° v are roughly equal, with the exception of a couple of anomalous points 
for surface UHMWl (the roughest surface) at low incidence angles. The backscat- 
tering coefficients seem somewhat noisier m the cross-polarized case. Because there 
was no test target that allowed easy and direct verification of measurement accuracy 
for cross-polarized returns, it is difficult to state probable measurement errors for 

“d * 1 - 

7.2 Comparison to Physical Optics Model 

The Kirchhoff or physical optics model is one of the most widely used surface 
scattering models [59]. The Kirchhoff model represents the scattered field in terms 
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Figure 7.2: Measured cross- polarized backscattering coefficients, a° h and <t£ v , of the 
three artificial surfaces: Surface A (UHMW1), Dj = 2.355; Surface B 
(UHMW3), D/ = 2.2275; and Surface C (UHMW2), Dj = 2.1. 


of the tangential field on the rough surface. It makes use of the tangent plane ap- 
proximation to the surface field, in which the surface field at a point is approximated 
by the field which would be present if the rough surface were replaced by a planar 
surface tangential to the surface at that point [56]. The tangential plane approxi- 
mation is valid if the radius of curvature at every point on the surface is large with 
respect to the radar wavelength. Use of the Kirchhoff model is usually limited to 
small incidence angles [14]. 

An additional simplification usually performed to make the physical optics so- 
lution more tractable is obtained by expanding the integrand of the scattered field 
integral into a series in terms of the surface slope and keeping only the lowest-order 
terms. This simplification is called the scalar approximation because it reduces to 
the sc alar formulation of Beckmann and Spizzichino [5]; it requires that the radius 
of curvature be larger than the radar wavelength and that the rms surface slope be 
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small. Commonly cited constraints are [56]: 


kl > 6 


R c > A 


m < 0.25 


where It is the wavenumber of the incident field, / is the surface correlation length. 
Rc is the average radius of curvature of the surface, and m is the rms slope. 

The scattering coefficient under the scalar approximation is given by [59] 

<4 = Cc + + C < 7 * 4 > 

where the subscripts p and q indicate the scattered and incident polarizations and 
the subscripts c, n, and s indicate the coherent, non-coherent, and slope terms. For 
backscatter at non-normal incidence, we omit the coherent term. 

The non-coherent term is given for the backscatter case by 

,0 = L^e— 


£ (4fc 2 cos* l^{-2k sin fl, 0) 


< 7 ° = 0 

pqn 


where R±o and R\\ 0 are the Fresnel reflection coefficients for horizontal and vertical 
polarization, respectively, and 

(7 - 7) 

For the power-law case, jRz(u,v) was calculated numerically from the modified 
roughness spectrum given by (5.8) using an inverse fast Fourier transform (FFT). 
The result was raised to the power n and then re- transformed to obtain ij^K 
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The slope term cr pp , was omitted for the power-law case because the non-analvtical 
correlation function prevented the direct evaluation of the slope terms. The a pqj is 
much smaller than in the Gaussian case, so it is assumed that the term is of 
little consequence. 

The physical optics results are compared with the measured <7° in Figure 7.3. 
The model values underestimate the measured backscatter by about 5 dB in the HH 
case, although they do follow the trend of the measured values out to about 45° or 
so. 

The W predictions start about 5 dB low at 15°, but rapidly grow worse due to the 
Brewster angle effect which the model has but which is not found in the measured 
data. The first-order physical optics formulation does not predict cross- polarized 
backscatter. 

The roughness spectra used in the physical optics power-law solution had a low- 
frequency cutoff of fu = (10A) -1 , where A is the radar wavelength. The low-frequency 
cutoff was imposed due to programming constraints in the calculation of the l£) 
terms. Letting fu vary from (5A) -1 to (25A) -1 resulted in a change in <r° of less than 
1 dB at 15° and still less at higher incidence angles. 

7.3 Comparison to Geometric Optics Model 

An alternative simplification to the Kirchhoff formulation involves the use of 
the stationary phase approximation, which means that scattering can occur only 
along directions for which there are specular points on the surface [59]. This high- 
frequency solution is valid when the average radius of curvature and the surface 
standard deviation are large compared to the radar wavelength. Commonly cited 
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Figure 7.3: Physical optics results for (a) HH and (b) W backscatter compared to 
the measured values for the three artificial surfaces. 
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criteria for this solution axe [56] 


kl > 6 


(7.8) 


R c > A 


ka > 


yTo 

cos 0, — cos 6 S 


(7.9) 


(7.10) 


where a is the standard deviation of surface height, 9, and 9, are the incident and 
scattered angles, and k, /, R c , and A are defined as before. 

The geometrical optics solution is given by Ulaby et al. [59]: 


, _ \mi 

w> 2 cos 4 0<7 2 |p"(O)| 


(7.11) 

(7.12) 


where R is the Fresnel reflection coefficient and p"( 0) is the second derivative of the 
normalized correlation function, p(t), 


P( T ) = 


Rzjr) 


(7.13) 


evaluated at zero. 

It is the p"(0) term that provides the main difficulty in applying (7.11) to the 
power-law case. Because the power-law correlation function has no analytical form, 
it must be evaluated numerically at discrete points. A discrete approximation was 
used to evaluate p"(0): 


_ _L _ El ZJ* \ 

- AV a a ) 

= £j0>, - 2* + Co) < 7 ' 14 > 

where the first three correlation values po, Pi, “d P* numerically calculated 
using the same bandlimited spectrum used for the physical optics case. For the three 




128 


7 



Figure 7.4: Geometrical optics results for backscattering compared to the measured 
values for the three artificial surfaces. 

artificial surfaces, A = 1.23 mm, and estimates of the rms slope, m = ^<7 2 |p"(0)|, 
were 0.31, 0.13, and 0.20. 

The geometrical optics results are compared with the measured backscatter in 
Figure 7.4. The model performs quite poorly for the power-law surfaces. 

7.4 Comparison to Small Perturbation Model 

The small perturbation model was first applied to electromagnetic scattering by 
Rice [49], who used the Rayleigh hypothesis to represent the scattered and trans- 
mitted fields near a rough surface in terms of waves traveling away from the surface 
only. These fields are expanded in a power series about some small parameter such 
as surface height or slope. The fields are then calculated to some low order, usually 
first order, by neglecting series terms of higher order [56]. Validity limits on the 
small perturbation model are 


ka < 0.3, 


(7.15) 
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m < 0.3, (7-16) 

where o. the standard deviation of surface elevation, is computed “from frequency 
components of the surface responsible for scattering at a given electromagnetic wave- 
length " [59]. The small perturbation model is usually used for larger incidence angles 
(0, > 20°) [14]. In a study of non-Gaussian surfaces by Chen et al. [15], the authors 
found that the small perturbation model has a wider range of validity when applied 
to a surface with a non-Gaussian (in their case, exponential) correlation function for 
large angles of incidence. 

The first-order small perturbation model solution for the backscatter case is given 
by 

= 81tV 3 cos' 4 9\a„\ 2 W(2ksm 0,0), (7.17) 

where 

Cthh = #i.(0)i 

sin 2 0 - Cr(l + sin 3 B) 

a w - (^- 1 )| £rCoe ^ + ( Cr _ ain *0)i/a]a’ 

&vh = = 0) 

and W(k x , k y ) is the normalized roughness spectrum, 

W{k z ,k y ) = ^-f [ p(u,u)exp[-j(**u + fc,»)]dudv (7.18) 

2* J—oo J — oo 

which is written in terms of the normalized surface correlation, p(u,v). The Fourier 
transform in (7.18) uses a different form than that used elsewhere in this dissertation, 
where the non-normalized correlation function, Rz( T *, T »)> i* defined as 

Rz(t x ,t v ) = r r° 5z(/,,/ f )exp[;2x(/,T, + / t T f )]d/ r d/ f . (7.19) 

«/— oo «/— OO 

The normalized spectrum W(Jfc», kg) with kg and kg in rad/m is therefore 

"sb 5 * (£■£)■ 


(7.20) 
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and we see that 

1 „ /Jtsin# 01 v 

W(2fcsinM) = 2 ^i Sz — '°J ’ ('* 21 ^ 

and therefore 

oj, = Sk* cos 4 6\a„\ 2 ^Sz ,0^ , (7.22) 

where the spectrum, Sz , is given by (5.8) or an equivalent expression. Note that the 
a 2 term drops out of the final SPM equation. 

The SPM results are shown in Figure 7.5 along with the measured backscatter. 
The SPM values track the measured data trends fairly well but have a negative 
offset of about 5 dB. In spite of this offset, the SPM results are the most useful of 
any of the models examined thus far. Because the SPM values are calculated based 
on a single frequency of the roughness spectrum, the offset cannot be attributed to 
an incorrect bandwidth choice. Because an error in the dielectric measurement was 
possible, different values of c, were tried. Changing c, from 2.05 to 2.3 raised <x° by 
about 1 dB. Errors in the measurement of more than 0.3 are considered extremely 
unlikely. Adding an imaginary component of c? = 2 raises a 2 by about 5 dB, but it 
is difficult to believe that the UHMW or the dielectric filler has a lossy part that is so 
radically different from tabulated values. In the absence of other explanations, the 
offset in the SPM predictions is viewed as a limitation of the model in the power-law 
case. 

7.5 Comparison to PWH model 

The last model examined is the Phased Wiener- Hermite (PWH) model for di- 
electric interfaces by Eftimiu [16], which is an extension of earlier work by Ef- 
timiu [17, 18, 19] and Eftimiu and Pan [20]. The model, which is relatively recent 
and which had no experimental verification, is based on an expansion of the surface 
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current in a series of orthogonal Wiener-Hermite functionals. Earlier uses of W iener- 
Hermite functional expansions in the study of turbulence led to their application to 
the problem of rough surface scattering by Nakayamaet al. [41, 42, 43] and Meecham 
and Lin [38], among others. The expansion is attractive because it represents the 
surface current in terms of functionals based on the rough surface random process, 
£(x,y), and because lower-order terms in the expansion contain reflections of all 
orders [38]. An overview of the PWH model is the subject of Appendix C. 

The PWH formulation is algebraically intensive in comparison with the other 
models in this chapter, perhaps because it relies on few simplifying assumptions. 
Aside from use of the Wiener-Hermite functional expansion (and the assumption 
that a first-order expansion is adequate), Eftimiu’s formulation uses the extinction 
theorem and the tangentiality of surface currents. The formulation is shown to reduce 
to the physical optics and small perturbation solutions in the appropriate limits. The 
PWH expression for o® is rather lengthy and involves a number of intermediate terms, 
and is therefore left in the Appendix. 

The PWH solutions for <r° of the three surface analogues are shown for HH . 
and W polarizations in Figure 7.6. The modeled which were calculated using 
the bandlimited roughness spectrum truncated at /j c = (10A) -1 , display a more- 
or-less linear dependence on incidence angle and predict higher backscattering for 
the rougher surfaces. However, the theoretical curves do not match the measured 
values very well, either in magnitude, difference between the rougher and smoother 
surfaces, or in angular trends. Recalculations in which a wider roughness spectrum 
(down to (25A) -1 ) was used resulted in rough and smooth a 0 values that were closer 
together in magnitude, but the values overestimated the measured values worse than 
in the narrower spectrum case. Values of a Q m were closer, differing from the measured 
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values by 6 dB or less. The trends of a ^ vs. 6 are closer than those of cr° K , but still 
not as close as those predicted by the small perturbation method. 

While there are no obvious flaws in the PWH approach, one can imagine sev- 
eral possible reasons why the model performs poorly in the power-law case. The 
algebraic complexity of the model is quite high, increasing the probability of math- 
ematical or programming errors. (The cited articles by Eftimiu contained numerous 
typographical errors; the author was supplied with a copy of Eftimiu's derivation 
notes to aid in reproducing the derivations.) The model requires the evaluation of a 
number of double integrals. These may be analytically reduced to one dimension in 
the Gaussian case but not in the power-law case, where the dynamic range (several 
orders of magnitude) may introduce problems in numerical precision. Finally, the 
power-law surface may simply be too complex to represent by a first-order Wiener- 
Hermite expansion, due to the presence of surface structure over a wide range of 
spatial scales. 
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(b) 

Figure 7.6: Phased Wiener- Hermite model results for (a) HH and (b) VV backseat - 
ter for the roughest and smoothest surfaces compared to the measured 
values. 
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VII Continuing Work 

Empirical laws or, to a limited extent, small perturbation theory can be used 
to relate (3 and c in equation (1) to scattering cross section. We are currently testing 
this approach by developing a classification algorithm from the empirical laws 
discussed in Section V. Using this algorithm, we will attempt to classify Andean 
Volcanics according to their roughness. This application is funded by a SIR-C 
project and was not a part of our Mars investigation. 
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